Hyunsoo Kim
  • Profile
  • Data Mining
  • R Basic
  • Data Visualization
  • Regression Analysis
  • Spatial Information Analysis
  • OpenData Analysis
  • Big Data Analysis Engineer

Regression Analysis

  • Profile
    • Profile
  • Study
    • Data Mining
    • R Basic
    • Data Visualization
    • Regression Analysis
    • Spatial Information Analysis
    • OpenData Analysis
    • Big Data Analysis Engineer

On this page

  • 0316 - 선형모형
    • 2장
      • 표2.5
      • 표2.6
      • 공분산(Covariance) - 식2.2
    • 0321 - 선형모형
      • 식 2.6
      • 그림 2.4
      • 표 2.3
      • 그림 2.2
      • 표 2.4
      • 그림 2.3
    • 단순선형회귀모형 2.4
      • 2.5 모수에 대한 추정
      • 식 2.14 & 2.15
      • 표2.7
      • 그림 2.5 - 산점도 + 회귀선
    • 0323
    • 0328
      • 2.7 신뢰구간
      • 2.8 예측
      • 2.9 적합성의 측정
      • 2.10 원점을 통과하는 회귀선
      • 2.11
    • 3장 다중선형회귀
      • 3.3 사례 감독자 직무수행능력 데이터
      • 3.4 모수 추정
      • 3.5 회귀계수에 대한 해석
      • 3.7 최소제곱추정량의 성질
      • 3.9 개별 회귀계수들에 대한 추론
      • 3.10 선형모형에서의 가설검정
      • 3.10.1 모든 회귀예수들이 0인가에 대한 검정
      • 3.10.2 회귀계수들의 부분집합이 0인가에 대한 검정
      • 3.10.3 회귀계수들의 통일성에 대한 검정
      • 3.10.4 제약조건하에서 회귀계수에 대한 추정과 검정
      • 3.11 예측
    • 부록 : 행렬을 이용한 회귀계수 추정
    • 4장 회귀진단: 모형위반의 검출
      • 4.3 다양한 유형의 잔차들
      • 4.5 모형을 적합깅 이전의 그래프
      • 4.5.1 일차원 그래프
      • 4.5.2 이차원 그래프
      • 4.7 선형성과 정규성 가정에 대한 검토
    • 0502 - 선형모형
      • 4.7 선형성과 정규성 가정에 대한 검토
      • Figure4.4
      • 4.8 지레점, 영향력, 특이값
      • 그림 [4.5]
      • 그림 [4.6] - 잔차의 인덱스플롯
      • 그림 [4.6] - 지레값의 인덱스플롯
      • 4.9 영향력의 측도
      • 4.9.1 Cook의 거리
      • 4.9.2 Welsch & Kuh의 측도
      • 4.9.3 Hadi의 영향력 측도
      • 4.10 잠재성 - 잔차플롯
      • 4.11 특이값에 대한 처리
      • 사례:스코틀랜드 언덕 경주 데이터
      • 4.14 로버스트 회귀 - outlier에 크게 영향을 받지 않음
    • 0509 - regression analysis
    • 5장 질적 예측 변수
    • 5.2 급료조사 데이터
      • 5.4 회귀방정식의 체계: 두 집단의 비교
      • 표 5.7 - 고용전 검사 프로그램 데이터
      • 그림 5.7 - 표준화잔차 대 검사점수 : 모형 1
      • 그림 5.8 - 표준화잔차 대 검사점수 : 모형 3
      • 그림 5.9 - 표준화잔차 대 검사점수 : 모형 1
      • 그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만
      • 그림 5.11 - 표준화잔차 대 검사점수 : 모형 1 . 백인만
      • 5.4.2 동일한 기울기와 다른 절편항을 가지는 모형
      • 5.4.3 동일한 절편항과 다른 기울기를 가지는 모형
    • 6장 변수변환 : Transformation of Varibables
      • 6.3 x-선 방사에 의한 박테리아 사망률
      • 그림 6.5
      • 6.3.1 선형모형의 부적절성
      • 그림 6.6
      • 6.3.2 선형성을 위한 로그변환
      • 그림 6.7
      • 6.4 분산안정화 변환
      • 그림 6.12
      • 6.5 이분산성의 검출
      • 6.6 이분산성의 제거
      • 6.7 가중최소제곱법
      • 6.8 데이터에 대한 로그변환
      • 9.2 통계적 추론에 미치는 효과
      • 그림 9.1
      • 9.3 예측에 미치는 효과
      • 그림 9.3 : 표준화잔차의 인덱스플롯
      • 그림 9.4 : 잔차플롯
      • 9.4 다중공선성의 탐색
      • 표 9.5 : 프랑스 경제 데이터
    • 10장 - 공선형 데이터의 처리
    • 11장
    • 0613 - regression analysis
      • 연습문제 11.5

Regression Analysis

Basic
Code
R
Author

Hyunsoo Kim

Published

March 14, 2022

0316 - 선형모형

예제를 통한 회귀분석

getwd() #"C:/Users/Hyunsoo Kim/Documents/lecture/regression_analysis"
[1] "C:/Users/Hyunsoo Kim/Desktop/senior_grade/blog/my-quarto-website"

2장

표2.5

data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

dim(data_2.5) #14 2
[1] 14  2
head(data_2.5)
  Minutes Units
1      23     1
2      29     2
3      49     3
4      64     4
5      74     4
6      87     5
X<-data_2.5$Units

Y<-data_2.5$Minutes

표2.6

df<-data.frame(

  #1:length(X),

  Y,

  X,

  Y-mean(Y),

  X-mean(X),

  (Y-mean(Y))^2,

  (X-mean(X))^2,

  (Y-mean(Y))*(X-mean(X))

)

df
     Y  X Y...mean.Y. X...mean.X. X.Y...mean.Y...2 X.X...mean.X...2
1   23  1 -74.2142857          -5     5.507760e+03               25
2   29  2 -68.2142857          -4     4.653189e+03               16
3   49  3 -48.2142857          -3     2.324617e+03                9
4   64  4 -33.2142857          -2     1.103189e+03                4
5   74  4 -23.2142857          -2     5.389031e+02                4
6   87  5 -10.2142857          -1     1.043316e+02                1
7   96  6  -1.2142857           0     1.474490e+00                0
8   97  6  -0.2142857           0     4.591837e-02                0
9  109  7  11.7857143           1     1.389031e+02                1
10 119  8  21.7857143           2     4.746173e+02                4
11 149  9  51.7857143           3     2.681760e+03                9
12 145  9  47.7857143           3     2.283474e+03                9
13 154 10  56.7857143           4     3.224617e+03               16
14 166 10  68.7857143           4     4.731474e+03               16
   X.Y...mean.Y......X...mean.X..
1                       371.07143
2                       272.85714
3                       144.64286
4                        66.42857
5                        46.42857
6                        10.21429
7                         0.00000
8                         0.00000
9                        11.78571
10                       43.57143
11                      155.35714
12                      143.35714
13                      227.14286
14                      275.14286

공분산(Covariance) - 식2.2

COV_XY<-sum((Y-mean(Y))*(X-mean(X))) / (length(X)-1) #136

### cov() 함수

cov(X,Y) #136
[1] 136
### 상관계수(correalationship)

### cor() 함수

cor(X,Y) #0.9936987 
[1] 0.9936987

0321 - 선형모형

data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

x<-data_2.5$Units

y<-data_2.5$Minutes

식 2.6

cor_xy<- COV_XY / (sd(x)*sd(y))

cor_xy
[1] 0.9936987
### cor() 함수

cor(x,y)
[1] 0.9936987
cor(y,x)
[1] 0.9936987
data_2.5
   Minutes Units
1       23     1
2       29     2
3       49     3
4       64     4
5       74     4
6       87     5
7       96     6
8       97     6
9      109     7
10     119     8
11     149     9
12     145     9
13     154    10
14     166    10
cor(data_2.5)
          Minutes     Units
Minutes 1.0000000 0.9936987
Units   0.9936987 1.0000000

그림 2.4

class(X)
[1] "numeric"
class(Y) #both numeric
[1] "numeric"
plot(X,Y, pch=19,xlab="Units",ylab="Minutes") 

표 2.3

data_2.3<-read.table("All_Data/p029a.txt",header=TRUE,sep="\t")

data_2.3
    Y  X
1   1 -7
2  14 -6
3  25 -5
4  34 -4
5  41 -3
6  46 -2
7  49 -1
8  50  0
9  49  1
10 46  2
11 41  3
12 34  4
13 25  5
14 14  6
15  1  7
X<-data_2.3$X

Y<-data_2.3$Y

그림 2.2

plot(X,Y)

cor(X,Y) # 0 완벽하게 2차함수의 형태도 0이 나옴(직선의 형태가 아닌것만)
[1] 0

표 2.4

data_2.4<-read.table("All_Data/p029b.txt",header=TRUE,sep="\t")

그림 2.3

plot(data_2.4$X1,data_2.4$Y1, pch=19); abline(3,0.5) #기울기 3 절편0.5인 선을 추가해라 

plot(data_2.4$X2,data_2.4$Y2, pch=19); abline(3,0.5)

plot(data_2.4$X3,data_2.4$Y3, pch=19); abline(3,0.5)

plot(data_2.4$X4,data_2.4$Y4, pch=19); abline(3,0.5)

m<-matrix(1:4,ncol=2,byrow=TRUE) #2행의 매트릭스 생성 

m
     [,1] [,2]
[1,]    1    2
[2,]    3    4
layout(m)

plot(Y1~X1, data=data_2.4,pch=19); abline(3,0.5)

#y~x  -> y=ax+b 이러한 형태를 가지는 모형식이라는 의미 

plot(Y2~X2, data=data_2.4,pch=19); abline(3,0.5)

plot(Y3~X3, data=data_2.4,pch=19); abline(3,0.5)

plot(Y4~X4, data=data_2.4,pch=19); abline(3,0.5)

layout(1) #변환을 다시 하지 않으면 설정한 매트릭스의 비율로 그래프가 그려짐 해제 필요 

# cor()

cor(data_2.4$X1,data_2.4$Y1) #0.8164205
[1] 0.8164205
cor(data_2.4$X2,data_2.4$Y2) #0.8162365
[1] 0.8162365
cor(data_2.4$X3,data_2.4$Y3) #0.8162867
[1] 0.8162867
cor(data_2.4$X4,data_2.4$Y4) #0.8165214
[1] 0.8165214
cor(data_2.4) #이렇게 한번에 할 수 있으나 가독성 떨어짐 
           Y1         X1         Y2         X2         Y3         X3         Y4
Y1  1.0000000  0.8164205  0.7500054  0.8164205  0.4687167  0.8164205 -0.4891162
X1  0.8164205  1.0000000  0.8162365  1.0000000  0.8162867  1.0000000 -0.3140467
Y2  0.7500054  0.8162365  1.0000000  0.8162365  0.5879193  0.8162365 -0.4780949
X2  0.8164205  1.0000000  0.8162365  1.0000000  0.8162867  1.0000000 -0.3140467
Y3  0.4687167  0.8162867  0.5879193  0.8162867  1.0000000  0.8162867 -0.1554718
X3  0.8164205  1.0000000  0.8162365  1.0000000  0.8162867  1.0000000 -0.3140467
Y4 -0.4891162 -0.3140467 -0.4780949 -0.3140467 -0.1554718 -0.3140467  1.0000000
X4 -0.5290927 -0.5000000 -0.7184365 -0.5000000 -0.3446610 -0.5000000  0.8165214
           X4
Y1 -0.5290927
X1 -0.5000000
Y2 -0.7184365
X2 -0.5000000
Y3 -0.3446610
X3 -0.5000000
Y4  0.8165214
X4  1.0000000

단순선형회귀모형 2.4

2.5 모수에 대한 추정

data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

x<-data_2.5$Units

y<-data_2.5$Minutes

식 2.14 & 2.15

sum((y-mean(y))*(x-mean(x))) #1768
[1] 1768
sum((x-mean(x))^2) #114
[1] 114
beta1_hat<-sum((y-mean(y))*(x-mean(x))) / sum((x-mean(x))^2)

beta1_hat #15.50877
[1] 15.50877
beta0_hat <- mean(y) - (beta1_hat*mean(x))

beta0_hat #4.161654
[1] 4.161654
### 최소제곱회귀 방정식

# Minutes = 4.161654 + 15.50877 * Units

plot(beta0_hat+beta1_hat*x,pch=19);

# 4개의 고장 난 부품을 수리하는데 걸리는 에측시간

4.161654 + 15.50877 * 4 #66.19673
[1] 66.19673
units<-4

beta0_hat + beta1_hat*units
[1] 66.19674
### 적합값(Fitted value)

y_hat<-beta0_hat + beta1_hat*(x)

### 최소 제곱 잔차(residual)

e<-y-y_hat

e #합이 0이라는 특징이 존재
 [1]  3.3295739 -6.1791980 -1.6879699 -2.1967419  7.8032581  5.2944862
 [7] -1.2142857 -0.2142857 -3.7230576 -9.2318296  5.2593985  1.2593985
[13] -5.2493734  6.7506266
sum(e) #1.278977e-13 0에 근사한 추지가 나옴
[1] 1.278977e-13

표2.7

df_2.7<-data.frame(

  x=x,

  y=y,

  y_hat,

  e

)

df_2.7
    x   y     y_hat          e
1   1  23  19.67043  3.3295739
2   2  29  35.17920 -6.1791980
3   3  49  50.68797 -1.6879699
4   4  64  66.19674 -2.1967419
5   4  74  66.19674  7.8032581
6   5  87  81.70551  5.2944862
7   6  96  97.21429 -1.2142857
8   6  97  97.21429 -0.2142857
9   7 109 112.72306 -3.7230576
10  8 119 128.23183 -9.2318296
11  9 149 143.74060  5.2593985
12  9 145 143.74060  1.2593985
13 10 154 159.24937 -5.2493734
14 10 166 159.24937  6.7506266
### lm() 함수 (linear model)

# Minutes = beta0 + beta1 * Units + epsilon

# 모형식 : y~x

lm(y~x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
      4.162       15.509  
res_lm<-lm(Minutes~Units,data=data_2.5)

res_lm

Call:
lm(formula = Minutes ~ Units, data = data_2.5)

Coefficients:
(Intercept)        Units  
      4.162       15.509  
# 리스트의 이름 

names(res_lm)
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        
# 회귀계수

res_lm$coefficients
(Intercept)       Units 
   4.161654   15.508772 
coef(res_lm)
(Intercept)       Units 
   4.161654   15.508772 
# 적합값

res_lm$fitted.values
        1         2         3         4         5         6         7         8 
 19.67043  35.17920  50.68797  66.19674  66.19674  81.70551  97.21429  97.21429 
        9        10        11        12        13        14 
112.72306 128.23183 143.74060 143.74060 159.24937 159.24937 
fitted(res_lm)
        1         2         3         4         5         6         7         8 
 19.67043  35.17920  50.68797  66.19674  66.19674  81.70551  97.21429  97.21429 
        9        10        11        12        13        14 
112.72306 128.23183 143.74060 143.74060 159.24937 159.24937 
# 최소제곱잔차

res_lm$residuals
         1          2          3          4          5          6          7 
 3.3295739 -6.1791980 -1.6879699 -2.1967419  7.8032581  5.2944862 -1.2142857 
         8          9         10         11         12         13         14 
-0.2142857 -3.7230576 -9.2318296  5.2593985  1.2593985 -5.2493734  6.7506266 
resid(res_lm)
         1          2          3          4          5          6          7 
 3.3295739 -6.1791980 -1.6879699 -2.1967419  7.8032581  5.2944862 -1.2142857 
         8          9         10         11         12         13         14 
-0.2142857 -3.7230576 -9.2318296  5.2593985  1.2593985 -5.2493734  6.7506266 
residuals(res_lm)
         1          2          3          4          5          6          7 
 3.3295739 -6.1791980 -1.6879699 -2.1967419  7.8032581  5.2944862 -1.2142857 
         8          9         10         11         12         13         14 
-0.2142857 -3.7230576 -9.2318296  5.2593985  1.2593985 -5.2493734  6.7506266 

그림 2.5 - 산점도 + 회귀선

plot(beta0_hat+beta1_hat*x,pch=19);

#abline(beta0_hat,beta1_hat)

abline(res_lm)

0323

data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

res_lm <- lm(Minutes ~ Units, data = data_2.5)

res_lm

Call:
lm(formula = Minutes ~ Units, data = data_2.5)

Coefficients:
(Intercept)        Units  
      4.162       15.509  
res_lm_summ<-summary(res_lm)

res_lm_summ #Pr(>|t|) - p-value

Call:
lm(formula = Minutes ~ Units, data = data_2.5)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2318 -3.3415 -0.7143  4.7769  7.8033 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    4.162      3.355    1.24    0.239    
Units         15.509      0.505   30.71 8.92e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.392 on 12 degrees of freedom
Multiple R-squared:  0.9874,    Adjusted R-squared:  0.9864 
F-statistic: 943.2 on 1 and 12 DF,  p-value: 8.916e-13
# unit은 시간에 영향을준다 약15.5분 만큼씩 

# coefficient에서 p-value에 대해서 알 수 있음 

# beta_0는 0이라고 보면되느냐? p-value가 0.05보다 크기에 

0328

2.7 신뢰구간

confint(res_lm) # beta_0,1의 95% 신뢰구간을 뽑아줌 
                2.5 %   97.5 %
(Intercept) -3.148482 11.47179
Units       14.408512 16.60903
?confint #level = 1-alpha
starting httpd help server ... done
confint(res_lm, level=0.90) # 90%의 신뢰구간
                 5 %     95 %
(Intercept) -1.81810 10.14141
Units       14.60875 16.40879

2.8 예측

# 4개의 고장난 부품을 수리하는 데에 걸리는 시간 예측

x<-4

4.161654 + 15.508772 *4
[1] 66.19674
res_lm$coefficients[1]+(res_lm$coefficients[2]*x)
(Intercept) 
   66.19674 
### predict()

df<-data.frame(Units=4) 

predict(res_lm,newdata=df) # res_lm을 만들때 사용한 데이터형식으로 만들어주어야함
       1 
66.19674 
res_lm_pred<-predict(res_lm,newdata=df,se.fit=TRUE)

### 예측값

res_lm_pred$fit
       1 
66.19674 
### 표준오차

res_lm_pred$se.fit # 평균반응에 대한 표준오차 
[1] 1.759688
### 예측한계

df<-data.frame(Units=4) #예제서는 4대기준

df<-data.frame(Units=1:10) #다른것도 보고 싶은경우 

res_lm_pred_int_p<-predict(res_lm,newdata=df,interval="prediction")

### 신뢰한계

df<-data.frame(Units=1:10) #다른것도 보고 싶은경우 

res_lm_pred_int_c<-predict(res_lm,newdata=df,interval="confidence") #둘의 차이를 보면 예측한계의 범위가 더큼 

### 예측한계 & 신뢰한계

# 신뢰한계는 평균에서 멀어지만 오차의범위가 커지고 평균에 다가갈수록 오차가 줄어듬

plot(Minutes~Units,data=data_2.5,pch=19)

abline(res_lm,col="red",lwd=2)

lines(1:10,res_lm_pred_int_p[,"lwr"],col="darkgreen")

lines(1:10,res_lm_pred_int_p[,"upr"],col="darkgreen")

lines(1:10,res_lm_pred_int_c[,"lwr"],col="blue")

lines(1:10,res_lm_pred_int_c[,"upr"],col="blue")

2.9 적합성의 측정

res_lm_summ<-summary(res_lm)

res_lm_summ #Multiple R-squared:0.9874 -> 반응변수의 전체변이중 98.94%가 예측변수에 의해 설명된다

Call:
lm(formula = Minutes ~ Units, data = data_2.5)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2318 -3.3415 -0.7143  4.7769  7.8033 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    4.162      3.355    1.24    0.239    
Units         15.509      0.505   30.71 8.92e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.392 on 12 degrees of freedom
Multiple R-squared:  0.9874,    Adjusted R-squared:  0.9864 
F-statistic: 943.2 on 1 and 12 DF,  p-value: 8.916e-13
# 만약 R-squared가 1이면 완벽한 선형의 관계 100%라는 것을 의미한다.

# R-squared는 변수가 들어갈수록 커지기에 adjust R-squared를 사용 추후 설명 

2.10 원점을 통과하는 회귀선

# Minutes = beta1 + Units + epsilon

res_lm_no<-lm(Minutes~Units-1,data=data_2.5)

res_lm_no

Call:
lm(formula = Minutes ~ Units - 1, data = data_2.5)

Coefficients:
Units  
16.07  
summary(res_lm_no)

Call:
lm(formula = Minutes ~ Units - 1, data = data_2.5)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.5955 -2.4733  0.4417  5.0243  9.7023 

Coefficients:
      Estimate Std. Error t value Pr(>|t|)    
Units  16.0744     0.2213   72.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.502 on 13 degrees of freedom
Multiple R-squared:  0.9975,    Adjusted R-squared:  0.9974 
F-statistic:  5274 on 1 and 13 DF,  p-value: < 2.2e-16
coef(summary(res_lm_no)) #rsquared=0.9975
      Estimate Std. Error  t value     Pr(>|t|)
Units 16.07443  0.2213341 72.62519 2.380325e-18

2.11

y<-rnorm(30)

t.test(y,mu=0)

    One Sample t-test

data:  y
t = -0.67175, df = 29, p-value = 0.5071
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.5134310  0.2595474
sample estimates:
 mean of x 
-0.1269418 
summary(lm(y~1))

Call:
lm(formula = y ~ 1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.96463 -0.67286  0.02797  0.69205  2.12145 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.1269     0.1890  -0.672    0.507

Residual standard error: 1.035 on 29 degrees of freedom

3장 다중선형회귀

3.3 사례 감독자 직무수행능력 데이터

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

dim(data_3.3)
[1] 30  7
class(data_3.3)
[1] "data.frame"
sapply(data_3.3,class) #all numeric
        Y        X1        X2        X3        X4        X5        X6 
"numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" 
summary(data_3.3) #모든변수가 numeric이면 분위수도 보여준다 
       Y               X1             X2              X3              X4       
 Min.   :40.00   Min.   :37.0   Min.   :30.00   Min.   :34.00   Min.   :43.00  
 1st Qu.:58.75   1st Qu.:58.5   1st Qu.:45.00   1st Qu.:47.00   1st Qu.:58.25  
 Median :65.50   Median :65.0   Median :51.50   Median :56.50   Median :63.50  
 Mean   :64.63   Mean   :66.6   Mean   :53.13   Mean   :56.37   Mean   :64.63  
 3rd Qu.:71.75   3rd Qu.:77.0   3rd Qu.:62.50   3rd Qu.:66.75   3rd Qu.:71.00  
 Max.   :85.00   Max.   :90.0   Max.   :83.00   Max.   :75.00   Max.   :88.00  
       X5              X6       
 Min.   :49.00   Min.   :25.00  
 1st Qu.:69.25   1st Qu.:35.00  
 Median :77.50   Median :41.00  
 Mean   :74.77   Mean   :42.93  
 3rd Qu.:80.00   3rd Qu.:47.75  
 Max.   :92.00   Max.   :72.00  
### 산점도 행렬

plot(data_3.3)

3.4 모수 추정

lm(Y~X1+X2+X3+X4+X5+X6,data=data_3.3)

Call:
lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data_3.3)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5  
   10.78708      0.61319     -0.07305      0.32033      0.08173      0.03838  
         X6  
   -0.21706  
lm(Y~.,data=data_3.3) # X1+X2+X3+X4+X5+X6쓰는 것이 아니라 .을 써서 모든 변수를 다써줌 

Call:
lm(formula = Y ~ ., data = data_3.3)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5  
   10.78708      0.61319     -0.07305      0.32033      0.08173      0.03838  
         X6  
   -0.21706  

3.5 회귀계수에 대한 해석

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

lm(Y~X1+X2,data=data_3.3)

Call:
lm(formula = Y ~ X1 + X2, data = data_3.3)

Coefficients:
(Intercept)           X1           X2  
   15.32762      0.78034     -0.05016  
# (Intercept)         X1           X2  

#  15.32762      0.78034     -0.05016

# 1) Y에서 X1 효과 제거

m1<-lm(Y~X1,data=data_3.3) # y prime

m1$residuals # x1이 설명하지 못한값 / x1의 효과를 제거한 값
           1            2            3            4            5            6 
 -9.86142016   0.32865220   3.80099328  -0.91673799   7.76411473 -12.87985944 
           7            8            9           10           11           12 
 -6.93517726   0.02794419  -4.25432454   6.59248165   9.62936020   7.34709147 
          13           14           15           16           17           18 
  7.83787183  -9.00893436   4.51872455  -1.29120309  -4.51815400   5.34709147 
          19           20           21           22           23           24 
 -2.19900672  -8.14368889   5.43928784   3.59248165 -11.18056744  -2.29688270 
          25           26           27           28           29           30 
  7.87475038  -6.48127545   7.02794419  -9.38907907   6.48184600   5.74567546 
# 2) X2에서 X1 효과 제거

m2<-lm(X2~X1,data=data_3.3) # x2 prime 

m2$residuals 
          1           2           3           4           5           6 
-15.1300345  -0.7994502  13.1223579  -6.2864182  -2.9818979   1.8178376 
          7           8           9          10          11          12 
-11.3385461  -7.4428019  10.9659742  -5.2603543   6.8439016  -2.7473223 
         13          14          15          16          17          18 
  6.2266138  21.4529422  -4.4688659 -15.1382816   1.4268783  15.2526777 
         19          20          21          22          23          24 
 -8.8776421  19.2787417  -6.4866827   1.7396457  -0.8255141   4.0524132 
         25          26          27          28          29          30 
 -4.6691304   7.5311341   0.5571981  -4.2082264   8.4268783 -22.0340258 
# 3) X1의 효과가 제거된 Y와 X2의 적합 - 원점을 지나는 회귀선

lm(m1$residuals~m2$residuals-1) # 원점을 지나면 -1를 하고 진행 // -3.25e-17

Call:
lm(formula = m1$residuals ~ m2$residuals - 1)

Coefficients:
m2$residuals  
    -0.05016  
# 다른 효과 없이(다른값이 고정) Y에 영향을 주는 순수한 X2의 값

# m2$residuals  : -0.05016  ==  X2 : -0.05016  

### 단위길이 척도화 - 잘사용하지않음

fn_scaling_len<-function(x){

  x0<-x-mean(x)

  x0/sqrt(sum(x0^2))

}

data_3.3_len<-sapply(data_3.3, fn_scaling_len)

data_3.3_len<-data.frame(data_3.3_len)

summary(data_3.3_len)
       Y                  X1                 X2                 X3           
 Min.   :-0.37579   Min.   :-0.41282   Min.   :-0.35109   Min.   :-0.353871  
 1st Qu.:-0.08975   1st Qu.:-0.11297   1st Qu.:-0.12344   1st Qu.:-0.148193  
 Median : 0.01322   Median :-0.02231   Median :-0.02479   Median : 0.002109  
 Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.000000  
 3rd Qu.: 0.10857   3rd Qu.: 0.14504   3rd Qu.: 0.14216   3rd Qu.: 0.164278  
 Max.   : 0.31070   Max.   : 0.32635   Max.   : 0.45328   Max.   : 0.294804  
       X4                 X5                 X6          
 Min.   :-0.38637   Min.   :-0.48356   Min.   :-0.32367  
 1st Qu.:-0.11401   1st Qu.:-0.10353   1st Qu.:-0.14318  
 Median :-0.02024   Median : 0.05130   Median :-0.03489  
 Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000  
 3rd Qu.: 0.11371   3rd Qu.: 0.09821   3rd Qu.: 0.08693  
 Max.   : 0.41733   Max.   : 0.32341   Max.   : 0.52461  
lm(Y~.,data=data_3.3_len)

Call:
lm(formula = Y ~ ., data = data_3.3_len)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5  
 -1.259e-16    6.707e-01   -7.343e-02    3.089e-01    6.981e-02    3.120e-02  
         X6  
 -1.835e-01  
### 표준화

# scale()

data_3.3_std<-scale(data_3.3)

#summary(data_3.3_std)

#sapply(data_3.3_std, sd, na.rm=T)

#class(data_3.3_std) #"matrix"

data_3.3_std<-data.frame(data_3.3_std)

#class(data_3.3_std) #"data.frame"

#sapply(data_3.3_std, sd, na.rm=T)

lm(Y~.,data=data_3.3_std) # beta게수 구하기 

Call:
lm(formula = Y ~ ., data = data_3.3_std)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5  
 -7.717e-16    6.707e-01   -7.343e-02    3.089e-01    6.981e-02    3.120e-02  
         X6  
 -1.835e-01  

3.7 최소제곱추정량의 성질

res_lm<-lm(Y~.,data=data_3.3)

res_lm

Call:
lm(formula = Y ~ ., data = data_3.3)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5  
   10.78708      0.61319     -0.07305      0.32033      0.08173      0.03838  
         X6  
   -0.21706  
summary(res_lm)

Call:
lm(formula = Y ~ ., data = data_3.3)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9418  -4.3555   0.3158   5.5425  11.5990 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.78708   11.58926   0.931 0.361634    
X1           0.61319    0.16098   3.809 0.000903 ***
X2          -0.07305    0.13572  -0.538 0.595594    
X3           0.32033    0.16852   1.901 0.069925 .  
X4           0.08173    0.22148   0.369 0.715480    
X5           0.03838    0.14700   0.261 0.796334    
X6          -0.21706    0.17821  -1.218 0.235577    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared:  0.7326,    Adjusted R-squared:  0.6628 
F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05
m1<-summary(lm(Y~X1+X2+X3+X4+X5+X6,data=data_3.3)) #Adjusted R-squared:  0.6628 

m2<-summary(lm(Y~X1+X2+X3+X4+X5,data=data_3.3)) #Adjusted R-squared:  0.6561 

# X6가 들어가는 것이 더 좋은 모델 

m1$adj.r.squared
[1] 0.662846
m2$adj.r.squared #summary에서 보다 더 정확하게 수치가 나옴 
[1] 0.6560539

3.9 개별 회귀계수들에 대한 추론

res_lm_summ<-summary(res_lm)

res_lm_summ #p-value의 존재는 무언가를 검정했다라는 반증

Call:
lm(formula = Y ~ ., data = data_3.3)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9418  -4.3555   0.3158   5.5425  11.5990 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.78708   11.58926   0.931 0.361634    
X1           0.61319    0.16098   3.809 0.000903 ***
X2          -0.07305    0.13572  -0.538 0.595594    
X3           0.32033    0.16852   1.901 0.069925 .  
X4           0.08173    0.22148   0.369 0.715480    
X5           0.03838    0.14700   0.261 0.796334    
X6          -0.21706    0.17821  -1.218 0.235577    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared:  0.7326,    Adjusted R-squared:  0.6628 
F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05
# p-value<0.05 H_1 귀무가설 채택 

# p-value>0.05 H_0 영가설 채택 // X1을 제외하고는 영가설 유의한 의미가 없음(Y에영향주는)

# 모두다 0이라는 가설을 가지고 분모 분자의 오차가 카이제곱을 따르고 거기서 나온 통계량

# F-분포 자유도는 분자 분모 두개를 가짐 //모아서 계산을 하기에 각각 계산하는것과 결과다름 

# 영가설-모든 회귀계수가 0이다.

# 대립가설-적어도 하나는 0이 아니다. p-value: 1.24e-05 <0.05 대립가설 채택 

# p-value가 0.05보다 작으면 대립가설 채택!!!!!! 기억해 

# 회귀계수에 대한 신뢰구간 - 95% 신뢰한계

confint(res_lm) #-13.18712881 ~ 34.7612816
                   2.5 %     97.5 %
(Intercept) -13.18712881 34.7612816
X1            0.28016866  0.9462066
X2           -0.35381806  0.2077178
X3           -0.02827872  0.6689430
X4           -0.37642935  0.5398936
X5           -0.26570179  0.3424647
X6           -0.58571106  0.1515977
#X1  0.28016866  0.9462066  사이에 0이 들어가있으면 영향을 준다라느걸 의미

#X2 -0.35381806  0.2077178  p-value없이도 알 수 있음 

#X5가 가장 영향이 적음 p-value가 가장 크기에(영향 효과의 크기를 비교할때)

#p-value가 작을 수록 영향을 많이 준다 beta값을 보는 것이 아닌 p-value를 보는 것 중요

#가장 의미있는 변수?->p-value가장 작은거 // 대립가설채택 Y에 영향을 가장

3.10 선형모형에서의 가설검정

3.10.1 모든 회귀예수들이 0인가에 대한 검정

# H_0: beta_1:beta_6=0

model_reduce<-lm(Y~1,data=data_3.3)

model_full<-lm(Y~.,data=data_3.3)

anova(model_reduce,model_full)
Analysis of Variance Table

Model 1: Y ~ 1
Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6
  Res.Df  RSS Df Sum of Sq      F   Pr(>F)    
1     29 4297                                 
2     23 1149  6      3148 10.502 1.24e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#대립가설 = 완전모형이 적절하다 / 1.24e-05 *** < 0.05 

#의미 있는 예측 변수가 한개 이상 존재한다 

summary(model_full) #summary에서 beta_1~beta_6까지 모두가 0이라는 가설로 진행을 이미함

Call:
lm(formula = Y ~ ., data = data_3.3)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9418  -4.3555   0.3158   5.5425  11.5990 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.78708   11.58926   0.931 0.361634    
X1           0.61319    0.16098   3.809 0.000903 ***
X2          -0.07305    0.13572  -0.538 0.595594    
X3           0.32033    0.16852   1.901 0.069925 .  
X4           0.08173    0.22148   0.369 0.715480    
X5           0.03838    0.14700   0.261 0.796334    
X6          -0.21706    0.17821  -1.218 0.235577    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared:  0.7326,    Adjusted R-squared:  0.6628 
F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05
#F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05

#예상// 가장의미있는변수? -> X1 이유?-> p-value 0.000903 로 가장 작기에 영향많이 줄것으로 예측 

3.10.2 회귀계수들의 부분집합이 0인가에 대한 검정

model_reduce<-lm(Y~X1+X3,data=data_3.3)

model_full<-lm(Y~.,data=data_3.3)

anova(model_reduce,model_full) #0.7158 > 0.05
Analysis of Variance Table

Model 1: Y ~ X1 + X3
Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     27 1254.7                           
2     23 1149.0  4    105.65 0.5287 0.7158
#영가설은 H_0: b_2=b_4=b_5=b_6=0 이라는 사실을 알 수 있다 

#b_1&b_3는 반응변수에 유의한 반응을 준다라는 것도 연계하여 알 수 있다 

3.10.3 회귀계수들의 통일성에 대한 검정

#해당 조건이 주어지고 만족할 때 beta_1=beta_3은 맞는가?

model_reduce<-lm(Y~I(X1-X3),data=data_3.3) #I를 씌우면 새로운 변수를 만든것과 동일

# X1-X3를 한 그자체를 분석하라는 의미//본래는 X1-X3 해서 새로운 변수를 만들어서 해야함 

model_full<-lm(Y~X1+X3,data=data_3.3)

anova(model_reduce,model_full) 
Analysis of Variance Table

Model 1: Y ~ I(X1 - X3)
Model 2: Y ~ X1 + X3
  Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
1     28 3846.7                                 
2     27 1254.6  1      2592 55.78 4.925e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#install.packages("car")

library(car)
Warning: package 'car' was built under R version 4.3.2
Loading required package: carData
Warning: package 'carData' was built under R version 4.3.2
model_full<-lm(Y~X1+X3,data=data_3.3)

car::linearHypothesis(model_full,c("X1=X3"))
Linear hypothesis test

Hypothesis:
X1 - X3 = 0

Model 1: restricted model
Model 2: Y ~ X1 + X3

  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     28 1424.6                              
2     27 1254.7  1    169.95 3.6572 0.06649 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.10.4 제약조건하에서 회귀계수에 대한 추정과 검정

# H_0: beta_1+beta_3=1 | beta_2=beta_3:beta_6=0

model_full<-lm(Y~X1+X3,data=data_3.3)

car::linearHypothesis(model_full,c("X1 + X3 = 1"))
Linear hypothesis test

Hypothesis:
X1  + X3 = 1

Model 1: restricted model
Model 2: Y ~ X1 + X3

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     28 1329.5                           
2     27 1254.7  1    74.898 1.6118 0.2151
# x1의 효과가 증가하면 x3의 효과는 감소한다 상대적인 관계 (반대로도 가능)

3.11 예측

model_full<-lm(Y~.,data_3.3)

# 예측값 - 적합값

model_full$fitted.values
       1        2        3        4        5        6        7        8 
51.11030 61.35277 69.93944 61.22684 74.45380 53.94185 67.14841 70.09701 
       9       10       11       12       13       14       15       16 
79.53099 59.19846 57.92572 55.40103 59.58168 70.21401 76.54933 84.54785 
      17       18       19       20       21       22       23       24 
76.15013 61.39736 68.01656 55.62014 42.60324 63.81902 63.66400 44.62475 
      25       26       27       28       29       30 
57.31710 67.84347 75.14036 56.04535 77.66053 76.87850 
# 예측한계(Prediction Limits)

predict(model_full,newdata = data_3.3,interval = "prediction")
        fit      lwr       upr
1  51.11030 34.16999  68.05060
2  61.35277 46.34536  76.36018
3  69.93944 53.94267  85.93622
4  61.22684 45.44586  77.00783
5  74.45380 59.17630  89.73129
6  53.94185 37.21813  70.66557
7  67.14841 51.64493  82.65189
8  70.09701 54.54384  85.65017
9  79.53099 62.71383  96.34814
10 59.19846 44.03506  74.36185
11 57.92572 42.00674  73.84470
12 55.40103 39.79333  71.00873
13 59.58168 43.39853  75.76483
14 70.21401 52.06636  88.36167
15 76.54933 60.79444  92.30422
16 84.54785 66.41374 102.68197
17 76.15013 59.99991  92.30036
18 61.39736 43.23384  79.56088
19 68.01656 52.44673  83.58639
20 55.62014 39.63744  71.60284
21 42.60324 26.35046  58.85603
22 63.81902 48.40145  79.23659
23 63.66400 48.56222  78.76578
24 44.62475 27.25435  61.99514
25 57.31710 41.29380  73.34041
26 67.84347 49.98605  85.70089
27 75.14036 59.31975  90.96097
28 56.04535 40.18723  71.90348
29 77.66053 61.97564  93.34541
30 76.87850 60.27441  93.48258
# 신뢰한계(Confidence limits)

predict(model_full,newdata = data_3.3,interval = "confidence")
        fit      lwr      upr
1  51.11030 42.55502 59.66557
2  61.35277 57.97029 64.73524
3  69.93944 63.44979 76.42909
4  61.22684 55.28897 67.16471
5  74.45380 70.02428 78.88332
6  53.94185 45.82386 62.05984
7  67.14841 61.99316 72.30367
8  70.09701 64.79421 75.39980
9  79.53099 71.22222 87.83975
10 59.19846 55.18008 63.21683
11 57.92572 51.63028 64.22116
12 55.40103 49.94035 60.86171
13 59.58168 52.64531 66.51805
14 70.21401 59.46431 80.96372
15 76.54933 70.68118 82.41748
16 84.54785 73.82102 95.27468
17 76.15013 69.29093 83.00933
18 61.39736 50.62090 72.17383
19 68.01656 62.66507 73.36805
20 55.62014 49.16527 62.07502
21 42.60324 35.50593 49.70055
22 63.81902 58.92819 68.70985
23 63.66400 59.88479 67.44321
24 44.62475 35.24662 54.00288
25 57.31710 50.76233 63.87187
26 67.84347 57.59134 78.09561
27 75.14036 69.09798 81.18275
28 56.04535 49.90540 62.18531
29 77.66053 71.98300 83.33806
30 76.87850 69.00992 84.74707

부록 : 행렬을 이용한 회귀계수 추정

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

Y<-data_3.3$Y

X<-data_3.3[,-1]

X<-cbind(1,X)

X<-as.matrix(X)

#beta_hat<-solve(t(X) %*% X) %*% t(X) %*% Y # %*%행렬 계산 

P<-solve(t(X) %*% X) %*% t(X)

beta_hat<- P %*% Y

lm(Y~.,data_3.3)

Call:
lm(formula = Y ~ ., data = data_3.3)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5  
   10.78708      0.61319     -0.07305      0.32033      0.08173      0.03838  
         X6  
   -0.21706  

4장 회귀진단: 모형위반의 검출

4.3 다양한 유형의 잔차들

# 표준화 잔차

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

res_lm<-lm(Y~.,data=data_3.3)

class(res_lm)
[1] "lm"
mode(res_lm)
[1] "list"
names(res_lm)
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        
res_lm$fitted.values
       1        2        3        4        5        6        7        8 
51.11030 61.35277 69.93944 61.22684 74.45380 53.94185 67.14841 70.09701 
       9       10       11       12       13       14       15       16 
79.53099 59.19846 57.92572 55.40103 59.58168 70.21401 76.54933 84.54785 
      17       18       19       20       21       22       23       24 
76.15013 61.39736 68.01656 55.62014 42.60324 63.81902 63.66400 44.62475 
      25       26       27       28       29       30 
57.31710 67.84347 75.14036 56.04535 77.66053 76.87850 
str(res_lm)
List of 12
 $ coefficients : Named num [1:7] 10.7871 0.6132 -0.0731 0.3203 0.0817 ...
  ..- attr(*, "names")= chr [1:7] "(Intercept)" "X1" "X2" "X3" ...
 $ residuals    : Named num [1:30] -8.11 1.647 1.061 -0.227 6.546 ...
  ..- attr(*, "names")= chr [1:30] "1" "2" "3" "4" ...
 $ effects      : Named num [1:30] -354.011 54.107 2.742 11.715 -0.971 ...
  ..- attr(*, "names")= chr [1:30] "(Intercept)" "X1" "X2" "X3" ...
 $ rank         : int 7
 $ fitted.values: Named num [1:30] 51.1 61.4 69.9 61.2 74.5 ...
  ..- attr(*, "names")= chr [1:30] "1" "2" "3" "4" ...
 $ assign       : int [1:7] 0 1 2 3 4 5 6
 $ qr           :List of 5
  ..$ qr   : num [1:30, 1:7] -5.477 0.183 0.183 0.183 0.183 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:7] "(Intercept)" "X1" "X2" "X3" ...
  .. ..- attr(*, "assign")= int [1:7] 0 1 2 3 4 5 6
  ..$ qraux: num [1:7] 1.18 1 1.29 1.1 1.07 ...
  ..$ pivot: int [1:7] 1 2 3 4 5 6 7
  ..$ tol  : num 1e-07
  ..$ rank : int 7
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 23
 $ xlevels      : Named list()
 $ call         : language lm(formula = Y ~ ., data = data_3.3)
 $ terms        :Classes 'terms', 'formula'  language Y ~ X1 + X2 + X3 + X4 + X5 + X6
  .. ..- attr(*, "variables")= language list(Y, X1, X2, X3, X4, X5, X6)
  .. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:7] "Y" "X1" "X2" "X3" ...
  .. .. .. ..$ : chr [1:6] "X1" "X2" "X3" "X4" ...
  .. ..- attr(*, "term.labels")= chr [1:6] "X1" "X2" "X3" "X4" ...
  .. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(Y, X1, X2, X3, X4, X5, X6)
  .. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
  .. .. ..- attr(*, "names")= chr [1:7] "Y" "X1" "X2" "X3" ...
 $ model        :'data.frame':  30 obs. of  7 variables:
  ..$ Y : num [1:30] 43 63 71 61 81 43 58 71 72 67 ...
  ..$ X1: num [1:30] 51 64 70 63 78 55 67 75 82 61 ...
  ..$ X2: num [1:30] 30 51 68 45 56 49 42 50 72 45 ...
  ..$ X3: num [1:30] 39 54 69 47 66 44 56 55 67 47 ...
  ..$ X4: num [1:30] 61 63 76 54 71 54 66 70 71 62 ...
  ..$ X5: num [1:30] 92 73 86 84 83 49 68 66 83 80 ...
  ..$ X6: num [1:30] 45 47 48 35 47 34 35 41 31 41 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language Y ~ X1 + X2 + X3 + X4 + X5 + X6
  .. .. ..- attr(*, "variables")= language list(Y, X1, X2, X3, X4, X5, X6)
  .. .. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:7] "Y" "X1" "X2" "X3" ...
  .. .. .. .. ..$ : chr [1:6] "X1" "X2" "X3" "X4" ...
  .. .. ..- attr(*, "term.labels")= chr [1:6] "X1" "X2" "X3" "X4" ...
  .. .. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(Y, X1, X2, X3, X4, X5, X6)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
  .. .. .. ..- attr(*, "names")= chr [1:7] "Y" "X1" "X2" "X3" ...
 - attr(*, "class")= chr "lm"
# 잔차 

res_lm$residuals
          1           2           3           4           5           6 
 -8.1102953   1.6472337   1.0605589  -0.2268416   6.5462010 -10.9418499 
          7           8           9          10          11          12 
 -9.1484140   0.9029929  -7.5309862   7.8015424   6.0742817  11.5989723 
         13          14          15          16          17          18 
  9.4183197  -2.2140147   0.4506705  -3.5478519  -2.1501319   3.6026355 
         19          20          21          22          23          24 
 -3.0165587  -5.6201442   7.3967582   0.1809831 -10.6639999  -4.6247464 
         25          26          27          28          29          30 
  5.6828983  -1.8434727   2.8596385  -8.0453540   7.3394730   5.1215016 
resid(res_lm) #실제값에서 예측된 값을 뺸값
          1           2           3           4           5           6 
 -8.1102953   1.6472337   1.0605589  -0.2268416   6.5462010 -10.9418499 
          7           8           9          10          11          12 
 -9.1484140   0.9029929  -7.5309862   7.8015424   6.0742817  11.5989723 
         13          14          15          16          17          18 
  9.4183197  -2.2140147   0.4506705  -3.5478519  -2.1501319   3.6026355 
         19          20          21          22          23          24 
 -3.0165587  -5.6201442   7.3967582   0.1809831 -10.6639999  -4.6247464 
         25          26          27          28          29          30 
  5.6828983  -1.8434727   2.8596385  -8.0453540   7.3394730   5.1215016 
### 내적 표준화잔차

rstandard(res_lm)
          1           2           3           4           5           6 
-1.41498026  0.23955370  0.16744867 -0.03512080  0.97184596 -1.86133876 
          7           8           9          10          11          12 
-1.38317210  0.13709194 -1.29490454  1.14799070  0.95218982  1.76906521 
         13          14          15          16          17          18 
 1.51371017 -0.46212316  0.06961486 -0.73868563 -0.34446368  0.75418016 
         19          20          21          22          23          24 
-0.45861365 -0.88618779  1.19699287  0.02717120 -1.56184734 -0.85286680 
         25          26          27          28          29          30 
 0.89948517 -0.36581416  0.44430497 -1.25422677  1.12683185  0.85971512 
### 외적 표준화잔차

MASS::studres(res_lm)
          1           2           3           4           5           6 
-1.44835328  0.23458097  0.16386794 -0.03434974  0.97062209 -1.97526518 
          7           8           9          10          11          12 
-1.41280382  0.13413337 -1.31529351  1.15637546  0.95017640  1.86145176 
         13          14          15          16          17          18 
 1.56019127 -0.45407837  0.06809185 -0.73117411 -0.33776450  0.74689589 
         19          20          21          22          23          24 
-0.45059801 -0.88189556  1.20894332  0.02657438 -1.61559196 -0.84763116 
         25          26          27          28          29          30 
 0.89560731 -0.35881868  0.43641573 -1.27088888  1.13380428  0.85466249 
redsid_df<-data.frame(

  Y=data_3.3$Y,

  Y_hat=res_lm$fitted.values,

  resid=resid(res_lm),

  rstandard=rstandard(res_lm),

  studres=MASS::studres(res_lm)

)

redsid_df
    Y    Y_hat       resid   rstandard     studres
1  43 51.11030  -8.1102953 -1.41498026 -1.44835328
2  63 61.35277   1.6472337  0.23955370  0.23458097
3  71 69.93944   1.0605589  0.16744867  0.16386794
4  61 61.22684  -0.2268416 -0.03512080 -0.03434974
5  81 74.45380   6.5462010  0.97184596  0.97062209
6  43 53.94185 -10.9418499 -1.86133876 -1.97526518
7  58 67.14841  -9.1484140 -1.38317210 -1.41280382
8  71 70.09701   0.9029929  0.13709194  0.13413337
9  72 79.53099  -7.5309862 -1.29490454 -1.31529351
10 67 59.19846   7.8015424  1.14799070  1.15637546
11 64 57.92572   6.0742817  0.95218982  0.95017640
12 67 55.40103  11.5989723  1.76906521  1.86145176
13 69 59.58168   9.4183197  1.51371017  1.56019127
14 68 70.21401  -2.2140147 -0.46212316 -0.45407837
15 77 76.54933   0.4506705  0.06961486  0.06809185
16 81 84.54785  -3.5478519 -0.73868563 -0.73117411
17 74 76.15013  -2.1501319 -0.34446368 -0.33776450
18 65 61.39736   3.6026355  0.75418016  0.74689589
19 65 68.01656  -3.0165587 -0.45861365 -0.45059801
20 50 55.62014  -5.6201442 -0.88618779 -0.88189556
21 50 42.60324   7.3967582  1.19699287  1.20894332
22 64 63.81902   0.1809831  0.02717120  0.02657438
23 53 63.66400 -10.6639999 -1.56184734 -1.61559196
24 40 44.62475  -4.6247464 -0.85286680 -0.84763116
25 63 57.31710   5.6828983  0.89948517  0.89560731
26 66 67.84347  -1.8434727 -0.36581416 -0.35881868
27 78 75.14036   2.8596385  0.44430497  0.43641573
28 48 56.04535  -8.0453540 -1.25422677 -1.27088888
29 85 77.66053   7.3394730  1.12683185  1.13380428
30 82 76.87850   5.1215016  0.85971512  0.85466249

4.5 모형을 적합깅 이전의 그래프

4.5.1 일차원 그래프

a<-rnorm(100,70,10) #연속형 데이터

# 히스토그램 

hist(a)

hist(a,breaks=5) #범위를 조절 막대의 5번 자름 

# 줄기 잎 그림 

stem(a)

  The decimal point is 1 digit(s) to the right of the |

  3 | 8
  4 | 
  5 | 01146678999999
  6 | 001122223444456666666777888899999
  7 | 000111222333333445555667777778889
  8 | 011122333334557789
  9 | 2
stem(round(a)) #줄기잎을 그릴때는 반올림을 하고 항상 진행 

  The decimal point is 1 digit(s) to the right of the |

  3 | 8
  4 | 
  5 | 01146678999999
  6 | 001122223444456666666777888899999
  7 | 000111222333333445555667777778889
  8 | 011122333334557789
  9 | 2
stem(round(a),scale=2) #scale을 2배로 늘려라 5기준으로 반으로 잘라서 

  The decimal point is 1 digit(s) to the right of the |

  3 | 8
  4 | 
  4 | 
  5 | 0114
  5 | 6678999999
  6 | 0011222234444
  6 | 56666666777888899999
  7 | 00011122233333344
  7 | 5555667777778889
  8 | 011122333334
  8 | 557789
  9 | 2
# 모든데이터를 볼 수있는 장점 데이터가 많으면 구림 

# 점플롯

idx<-rep(1,length(a)) #a의 갯수에 맞춰서 1를 반복 

plot(idx,a)

plot(jitter(idx),a,xlim=c(0.5,1.5))

# 상자그림

boxplot(a) #사분위수에 대해서 알 수 있음 

# 상자그림 + 점플롯

boxplot(a)

points(jitter(idx),a)

4.5.2 이차원 그래프

data_4.1<-read.table("All_Data/p103.txt",header=T,sep="\t")

data_4.1
       Y   X1   X2
1  12.37 2.23 9.66
2  12.66 2.57 8.94
3  12.00 3.87 4.40
4  11.93 3.10 6.64
5  11.06 3.39 4.91
6  13.03 2.83 8.52
7  13.13 3.02 8.04
8  11.44 2.14 9.05
9  12.86 3.04 7.71
10 10.84 3.26 5.11
11 11.20 3.39 5.05
12 11.56 2.35 8.51
13 10.83 2.76 6.59
14 12.63 3.90 4.90
15 12.46 3.16 6.96
class(data_4.1) #data.frame
[1] "data.frame"
# 산점도 행렬

plot(data_4.1)

cor(data_4.1) #상관계수
             Y           X1         X2
Y  1.000000000  0.002497966  0.4340688
X1 0.002497966  1.000000000 -0.8997765
X2 0.434068758 -0.899776481  1.0000000
pairs(data_4.1)

# Correlation panel

panel.cor<-function(x,y){

  par(usr=c(0,1,0,1))

  r<-round(cor(x,y),digits = 3)

  text(0.5,0.5,r,cex=1.5)

}

pairs(data_4.1,lower.panel = panel.cor)

# 회전도표, 동적 그래프(3차원)

#install.packages("rgl")

library(rgl)
Warning: package 'rgl' was built under R version 4.3.2
plot3d(x=data_4.1$X1,y=data_4.1$X2,z=data_4.1$Y) 

4.7 선형성과 정규성 가정에 대한 검토

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

res_lm<-lm(Y~X1+X3,data=data_3.3)

layout(matrix(1:4,nrow=2,byrow=T))

plot(res_lm) #회귀진단 플랏이 나옴 

layout(1)

# 1. 표준화잔차의 정규확률분포

plot(res_lm,2)

# 2. 표준화잔차 대 각 예측변수들의 산점도

plot(res_lm,3) #이런경우에는 데이터를 추가한다 좌측하단이 없어서 

# 3. 표준화잔차 대 적합값의 플롯

# 4. 표준화잔차의 인덱스 플롯

0502 - 선형모형

library(dplyr)

Attaching package: 'dplyr'
The following object is masked from 'package:car':

    recode
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

4.7 선형성과 정규성 가정에 대한 검토

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

res_lm<-lm(Y~X1+X3,data=data_3.3) #두개의 변수만 의미있다고 가정하고 진행 

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(res_lm) #회귀진단 플랏이 나옴 / 총6개임 

layout(1) #다시 한개의 플랏만 그리도록 

# 1. 표준화잔차의 정규확률분포

plot(res_lm,2) #QQ-plot y=x 기울기의 직선위에 점들이 있어야 한다 눈대중으로 

# 2. 표준화잔차 대 각 예측변수들의 산점도

plot(res_lm,3) #이런경우에는 데이터를 추가한다 좌측하단이 없어서 

#랜덤하게 데이터가 흩어져 있어야 한다

# 내적 표준화잔차

plot(data_3.3$X1,rstandard(res_lm))

plot(data_3.3$X3,rstandard(res_lm)) #각각의 잔차들이 랜덤하게 잘 퍼져야함 

# 3. 표준화잔차 대 적합값의 플롯

plot(res_lm,1) #잔차와 적합값은 상관성이 없어야하며 랜덤하게 퍼져야함 

# 4. 표준화잔차의 인덱스 플롯 

plot(res_lm,5)

data_2.4<-read.table("All_Data/p029b.txt",header=T,sep="\t")

m<-matrix(1:4,ncol=2,byrow=TRUE)  

layout(m)

plot(Y1~X1, data=data_2.4,pch=19); abline(3,0.5)

plot(Y2~X2, data=data_2.4,pch=19); abline(3,0.5)

plot(Y3~X3, data=data_2.4,pch=19); abline(3,0.5)

plot(Y4~X4, data=data_2.4,pch=19); abline(3,0.5)

layout(1)

Figure4.4

res_lm<-lm(Y1~X1,data=data_2.4) 

#1번플랏은 적당히 잘퍼짐,2번플랏은 어느정도 선형성 있음(데이터적어서그런거임)

res_lm<-lm(Y2~X2,data=data_2.4)

res_lm<-lm(Y3~X3,data=data_2.4)

res_lm<-lm(Y4~X4,data=data_2.4)

res_lm

Call:
lm(formula = Y4 ~ X4, data = data_2.4)

Coefficients:
(Intercept)           X4  
     3.0017       0.4999  
layout(matrix(1:4,nrow=2,byrow=T))

plot(res_lm) 
Warning: not plotting observations with leverage one:
  8

layout(1)

4.8 지레점, 영향력, 특이값

# 사례: 뉴욕 강 데이터

# agr-농업, forest-숲, rsdntial-주거, comlndl-산업, nitrogen-질소

data_1.9<-read.table("All_Data/p010.txt",header=T,sep="\t")

head(data_1.9)
       River Agr Forest Rsdntial ComIndl Nitrogen
1      Olean  26     63      1.2    0.29     1.10
2  Cassadaga  29     57      0.7    0.09     1.01
3      Oatka  54     26      1.8    0.58     1.90
4  Neversink   2     84      1.9    1.98     1.00
5 Hackensack   3     27     29.4    3.11     1.99
6  Wappinger  19     61      3.4    0.56     1.42
plot(data_1.9[-1],pch=19) #river라는 첫번째 컬럼을 제외하고 진행 

res_1<-lm(Nitrogen~.,data=data_1.9[-1]) #모든데이터 사용

res_2<-lm(Nitrogen~.,data=data_1.9[-4,-1]) #4번쨰 데이터 제외

res_3<-lm(Nitrogen~.,data=data_1.9[-5,-1]) #5번쨰 데이터 제외 

#회귀계수

data.frame(all=coef(res_1),

           rm4=coef(res_2),

           rm5=coef(res_3))
                     all          rm4          rm5
(Intercept)  1.722213529  1.099471134  1.626014115
Agr          0.005809126  0.010136685  0.002352222
Forest      -0.012967887 -0.007589231 -0.012760349
Rsdntial    -0.007226768 -0.123792917  0.181160986
ComIndl      0.305027765  1.528956204  0.075617570
#p-value

data.frame(all=coef(summary(res_1))[,4],

           rm4=coef(summary(res_2))[,4],

           rm5=coef(summary(res_3))[,4])
                   all          rm4        rm5
(Intercept) 0.18316946 0.2477883387 0.05619948
Agr         0.70462624 0.3717054741 0.80880737
Forest      0.36667966 0.4700975391 0.16975563
Rsdntial    0.83372002 0.0071342930 0.00112280
ComIndl     0.08230952 0.0005512227 0.51774981
#4,5번째 데이터를 각각 뺴고 진행을 해보니 영향을 끼치는 값임을 알수있고

#5번째는 주거 관련해서는 부호를 바꿀 정도로 강력하다 

# 단순선형회구모형

res<-lm(Nitrogen~ComIndl,data=data_1.9)

그림 [4.5]

plot(Nitrogen~ComIndl,data=data_1.9,pch=19)

abline(res)

#leverage values 지레값

p_ii<-hatvalues(res)

hiegh_leverage<-ifelse(p_ii>2*2/30,data_1.9$River,"")

hiegh_leverage #높은 지레값을 가지고 있는 강의 이름을 표시(이는 보기에 편하기 위해서함)
           1            2            3            4            5            6 
          ""           ""           ""  "Neversink" "Hackensack"           "" 
           7            8            9           10           11           12 
          ""           ""           ""           ""           ""           "" 
          13           14           15           16           17           18 
          ""           ""           ""           ""           ""           "" 
          19           20 
          ""           "" 
text(data_1.9$ComIndl,data_1.9$Nitrogen-0.1,hiegh_leverage)

그림 [4.6] - 잔차의 인덱스플롯

plot(rstandard(res),pch=19) #2또는 3시그마를 넘으면 특이값이라 함 

그림 [4.6] - 지레값의 인덱스플롯

plot(p_ii,pch=19) #평균의 2배를 기준으로 비교함 

abline(h=2*2/30,col="red") #이보다 높은것이 지레값이 높은것 높은 영향력을 가진것 

# 단순선형회귀모형

res<-lm(Nitrogen~ComIndl,data=data_1.9)

plot(Nitrogen~ComIndl,data=data_1.9,pch=19)

abline(res)

res<-lm(Nitrogen~ComIndl,data=data_1.9[-4,-1])

plot(Nitrogen~ComIndl,data=data_1.9[-4,-1],pch=19)

abline(res) #4번째 제외

res<-lm(Nitrogen~ComIndl,data=data_1.9[-5,-1])

plot(Nitrogen~ComIndl,data=data_1.9[-5,-1],pch=19)

abline(res) #5번쨰 제외

#이처럼 4,5번을 빼고 진행을 하면 조금더 잘 나타냄

res<-lm(Nitrogen~ComIndl,data=data_1.9[c(-4,-5),-1])

plot(Nitrogen~ComIndl,data=data_1.9[-5,-1],pch=19)

abline(res) #4,5번째 제외 

4.9 영향력의 측도

4.9.1 Cook의 거리

res<-lm(Nitrogen~ComIndl,data=data_1.9)

plot(res,4)

#install.packages("olsrr")

library(olsrr)
Warning: package 'olsrr' was built under R version 4.3.2

Attaching package: 'olsrr'
The following object is masked from 'package:datasets':

    rivers

ols_plot_cooksd_chart(res)

4.9.2 Welsch & Kuh의 측도

olsrr::ols_plot_dffits(res)

4.9.3 Hadi의 영향력 측도

olsrr::ols_plot_hadi(res)

# Residual & Leverage & Cook's distance

plot(res,5) #영향력 관측치를 보기 위한 플랏 

4.10 잠재성 - 잔차플롯

olsrr::ols_plot_resid_pot(res) #2.0에 있는 값외에도 x축 0.2이상의 것들도 특이값으로 

#지레값이 커도 영향력이 없는 애들은 신경 안써도 되나 영향력이 큰애들을 탐색해 보아야 함

4.11 특이값에 대한 처리

#특이값이 큰경우 그것이 TRUE데이터면 오류가 있는 경우 수정을 하거나 가중치를 변화를

#하거나 데이터를 수정을 시켜주거나 다시 실험을 하는 등 여러가지 방법을 사용하여.. 

## 첨가변수플롯(added-variable plot), 편회귀플롯(partial regression plot)

res<-lm(Nitrogen~ComIndl,data=data_1.9)

car::avPlots(res,pch=19)

res<-lm(Nitrogen~.,data=data_1.9[-1])

car::avPlots(res,pch=19)

#beta별로 각각의 어떤 변수가 영향력을 많이 주는지 알게하는 함수 

사례:스코틀랜드 언덕 경주 데이터

data_4.5<-read.table("All_Data/p120.txt",header=T,sep="\t")

dim(data_4.5)
[1] 35  4
names(data_4.5)
[1] "Hill.Race" "Time"      "Distance"  "Climb"    
head(data_4.5)
                    Hill.Race Time Distance Climb
1 Greenmantle New Year Dash    965      2.5   650
2                  Carnethy   2901      6.0  2500
3              Craig Dunain   2019      6.0   900
4                   Ben Rha   2736      7.5   800
5                Ben Lomond   3736      8.0  3070
6                  Goatfell   4393      8.0  2866
# 회전도표 

library(rgl)

with(data_4.5,plot3d(x=Distance,y=Climb,z=Time))

res<-lm(Time~Distance+Climb,data=data_4.5)

summary(res) #beta_0가 마이너스여도 신경안씀 관심있는 것은 회귀계수임 

Call:
lm(formula = Time ~ Distance + Climb, data = data_4.5)

Residuals:
   Min     1Q Median     3Q    Max 
-973.0 -427.7  -71.2  142.2 3907.3 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -539.4829   258.1607  -2.090   0.0447 *  
Distance     373.0727    36.0684  10.343 9.86e-12 ***
Climb          0.6629     0.1231   5.387 6.44e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 880.5 on 32 degrees of freedom
Multiple R-squared:  0.9191,    Adjusted R-squared:  0.914 
F-statistic: 181.7 on 2 and 32 DF,  p-value: < 2.2e-16
#Time=-539.4829+373.0727*Distance+0.6629*Climb 
### 첨가변수플롯(added-variable plot), 편회귀플롯(partial regression plot)

car::avPlots(res,pch=19)

### 성분잔차플롯(component plus residual plots), 편자차플롯(partial residual plot)

car::crPlots(res,id=T,pch=19) #첨가변수보다 성분잔차를 더 많이 사용 / 비선형여부를 확인

# 점선에 비해서 분홍선이 크게 떨어져 있지 않아 선형적인 추세를 가지고 있다고 추정이 가능 

### 잠재성-잔차플롯

olsrr::ols_plot_resid_pot(res)

### Hadi의 영향력 측도

olsrr::ols_plot_hadi(res)

### Cook의 거리

olsrr::ols_plot_cooksd_chart(res) #전체적은 플롯을 보면서 이상치에 대한 확인을 함 

# 이러한 값들의 제외 여부는 연구자가 선택해서 진행을 함 

# outlier에 대한 처리를 어떻게 했다고 말을 해야함 

# 보고서는 옆에 사람이 보고 쉽게 따라할 수 있을 정도로 

# 어떤 속성으로 어떻게 진행을 했다라는 것을 중시 코드 보단 결과를 보여줘라 

4.14 로버스트 회귀 - outlier에 크게 영향을 받지 않음

library(MASS)

Attaching package: 'MASS'
The following object is masked from 'package:olsrr':

    cement
The following object is masked from 'package:dplyr':

    select
res<-lm(Time~Distance+Climb,data=data_4.5)

res

Call:
lm(formula = Time ~ Distance + Climb, data = data_4.5)

Coefficients:
(Intercept)     Distance        Climb  
  -539.4829     373.0727       0.6629  
res_rlm<-MASS::rlm(Time~Distance+Climb,data=data_4.5)

res_rlm
Call:
rlm(formula = Time ~ Distance + Climb, data = data_4.5)
Converged in 10 iterations

Coefficients:
 (Intercept)     Distance        Climb 
-576.3836570  393.0374614    0.4977894 

Degrees of freedom: 35 total; 32 residual
Scale estimate: 313 
summary(res)

Call:
lm(formula = Time ~ Distance + Climb, data = data_4.5)

Residuals:
   Min     1Q Median     3Q    Max 
-973.0 -427.7  -71.2  142.2 3907.3 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -539.4829   258.1607  -2.090   0.0447 *  
Distance     373.0727    36.0684  10.343 9.86e-12 ***
Climb          0.6629     0.1231   5.387 6.44e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 880.5 on 32 degrees of freedom
Multiple R-squared:  0.9191,    Adjusted R-squared:  0.914 
F-statistic: 181.7 on 2 and 32 DF,  p-value: < 2.2e-16
summary(res_rlm)

Call: rlm(formula = Time ~ Distance + Climb, data = data_4.5)
Residuals:
     Min       1Q   Median       3Q      Max 
-645.074 -197.082   -2.035  212.266 3942.045 

Coefficients:
            Value     Std. Error t value  
(Intercept) -576.3837  105.2774    -5.4749
Distance     393.0375   14.7086    26.7216
Climb          0.4978    0.0502     9.9200

Residual standard error: 312.6 on 32 degrees of freedom
#install.packages("robustbase")

library(robustbase)
Warning: package 'robustbase' was built under R version 4.3.2
res_lmrob<-lmrob(Time~Distance+Climb,data=data_4.5)

res_lmrob

Call:
lmrob(formula = Time ~ Distance + Climb, data = data_4.5)
 \--> method = "MM"
Coefficients:
(Intercept)     Distance        Climb  
  -487.3793     398.2784       0.3901  
summary(res_lmrob)

Call:
lmrob(formula = Time ~ Distance + Climb, data = data_4.5)
 \--> method = "MM"
Residuals:
    Min      1Q  Median      3Q     Max 
-650.01 -160.60   23.37  216.11 3875.00 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -487.37929   86.33384  -5.645 3.04e-06 ***
Distance     398.27843    5.98522  66.544  < 2e-16 ***
Climb          0.39013    0.04165   9.368 1.09e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 290.8 
Multiple R-squared:  0.9868,    Adjusted R-squared:  0.986 
Convergence in 9 IRWLS iterations

Robustness weights: 
 3 observations c(7,18,33) are outliers with |weight| = 0 ( < 0.0029); 
 2 weights are ~= 1. The remaining 30 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.5965  0.9098  0.9623  0.9214  0.9875  0.9990 
Algorithmic parameters: 
       tuning.chi                bb        tuning.psi        refine.tol 
        1.548e+00         5.000e-01         4.685e+00         1.000e-07 
          rel.tol         scale.tol         solve.tol          zero.tol 
        1.000e-07         1.000e-10         1.000e-07         1.000e-10 
      eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
        2.857e-03         1.364e-08         5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
           200              0           1000              0           2000 
                  psi           subsampling                   cov 
           "bisquare"         "nonsingular"         ".vcov.avar1" 
compute.outlier.stats 
                 "SM" 
seed : int(0) 
#standard error가 res > res_rlm > res_lmrob 순으로 되어 있음 

# 4장 연습문제 해보기

0509 - regression analysis

library(dplyr)

5장 질적 예측 변수

#install.packages("fastDummies")

library(fastDummies)
Warning: package 'fastDummies' was built under R version 4.3.2
Thank you for using fastDummies!
To acknowledge our work, please cite the package:
Kaplan, J. & Schlegel, B. (2023). fastDummies: Fast Creation of Dummy (Binary) Columns and Rows from Categorical Variables. Version 1.7.1. URL: https://github.com/jacobkap/fastDummies, https://jacobkap.github.io/fastDummies/.

5.2 급료조사 데이터

data_5.1<-read.table("All_Data/p130.txt",header=T,sep="\t")

# S:급료 X:경력 E:교육수준 M:관리(형태) / E,M은 범주형 변수

names(data_5.1)
[1] "S" "X" "E" "M"
# 범주형 질적 변수를 수치형으로 변형시켜서 예측하는데 사용한것이 질적 예측변수이다 

# E_1,E_2,E_3이런식으로 나누어서 0,1로 분류를 한다 (이것이 가변수)

# 더미변수를 만들 경우에는 역행렬의 조건에 의해서 -1개의 변수만 만들면 된다 

# 이는 공산성의 문제또한 있기에 이를 위해서 -1를 한것임 

### 자료형 변경 : 정수 -> 범주

data_5.1$E<-as.factor(data_5.1$E)

data_5.1$M<-as.factor(data_5.1$M)

head(data_5.1)
      S X E M
1 13876 1 1 1
2 11608 1 3 0
3 18701 1 3 1
4 11283 1 2 0
5 11767 1 3 0
6 20872 2 2 1
data_5.1$E #Levels: 1 2 3이라는 것이 생김 문자로 처리한다는 의미 
 [1] 1 3 3 2 3 2 2 1 3 2 1 2 3 1 3 3 2 2 3 1 1 3 2 2 1 2 1 3 1 1 2 3 2 2 1 2 3 1
[39] 2 2 3 2 2 1 2 1
Levels: 1 2 3
### 가변수 생성

data_5.1$E<-factor(as.character(data_5.1$E),levels = c("3","1","2")) 

#3번을 베이스 카테고리로 쓰기위한 설정 / 설정을 안하면 베이스는 e_1이 된다

data_5.1$M<-factor(as.character(data_5.1$M),levels = c("0","1")) 

data_dummy<-dummy_cols(data_5.1,

                       select_columns = c("E","M"),

                       remove_first_dummy = T,

                       remove_selected_columns = T) #첫번째 생성되는 더미 변수를 제거

data_dummy #가변수의 더미는 n-1개를 하는 것이 역행렬을 위한 것이기에 지워준다 
       S  X E_1 E_2 M_1
1  13876  1   1   0   1
2  11608  1   0   0   0
3  18701  1   0   0   1
4  11283  1   0   1   0
5  11767  1   0   0   0
6  20872  2   0   1   1
7  11772  2   0   1   0
8  10535  2   1   0   0
9  12195  2   0   0   0
10 12313  3   0   1   0
11 14975  3   1   0   1
12 21371  3   0   1   1
13 19800  3   0   0   1
14 11417  4   1   0   0
15 20263  4   0   0   1
16 13231  4   0   0   0
17 12884  4   0   1   0
18 13245  5   0   1   0
19 13677  5   0   0   0
20 15965  5   1   0   1
21 12336  6   1   0   0
22 21352  6   0   0   1
23 13839  6   0   1   0
24 22884  6   0   1   1
25 16978  7   1   0   1
26 14803  8   0   1   0
27 17404  8   1   0   1
28 22184  8   0   0   1
29 13548  8   1   0   0
30 14467 10   1   0   0
31 15942 10   0   1   0
32 23174 10   0   0   1
33 23780 10   0   1   1
34 25410 11   0   1   1
35 14861 11   1   0   0
36 16882 12   0   1   0
37 24170 12   0   0   1
38 15990 13   1   0   0
39 26330 13   0   1   1
40 17949 14   0   1   0
41 25685 15   0   0   1
42 27837 16   0   1   1
43 18838 16   0   1   0
44 17483 16   1   0   0
45 19207 17   0   1   0
46 19346 20   1   0   0
# 더미를 만든 이후 더미의 모체 변수인 E M을 지워주어야 한다 

### 회귀분석(1) - 가변수

res<-lm(S~.,data = data_dummy)

res 

Call:
lm(formula = S ~ ., data = data_dummy)

Coefficients:
(Intercept)            X          E_1          E_2          M_1  
    11031.8        546.2      -2996.2        147.8       6883.5  
### 회귀분석(2) - lm()

res_1<-lm(S~.,data=data_5.1)

res_1 

Call:
lm(formula = S ~ ., data = data_5.1)

Coefficients:
(Intercept)            X           E1           E2           M1  
    11031.8        546.2      -2996.2        147.8       6883.5  
#더미변수를 많이 쓰기에 factor로 바꾸어주고 분석하면 알아서 더미변수를 만들어서 진행함

summary(res_1)

Call:
lm(formula = S ~ ., data = data_5.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1884.60  -653.60    22.23   844.85  1716.47 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11031.81     383.22  28.787  < 2e-16 ***
X             546.18      30.52  17.896  < 2e-16 ***
E1          -2996.21     411.75  -7.277 6.72e-09 ***
E2            147.82     387.66   0.381    0.705    
M1           6883.53     313.92  21.928  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1027 on 41 degrees of freedom
Multiple R-squared:  0.9568,    Adjusted R-squared:  0.9525 
F-statistic: 226.8 on 4 and 41 DF,  p-value: < 2.2e-16
#E_3와 M_0는 Intercept에 포함이 되어있음 그렇기에 베이스 카테고리라고 한다 

#E가 3이고 M이 0이면 대학원이상 일반관리 직급 -> Intercept(Beta_0) + X*Beta_1 = Y

### 표준화잔차 대 경력연수

plot(data_5.1$X, rstandard(res_1),pch=19,xlab="범주",ylab="잔차")

# 0을 중심으로 잘 퍼져있는가를 확인해야함 

#### 표준화잔차 대 교육수준 - 관리 조합

EM<-paste0(data_5.1$E,data_5.1$M)

plot(EM,rstandard(res_1),pch=19,xlabl="범주",ylab="잔차")
Warning in plot.window(...): "xlabl" is not a graphical parameter
Warning in plot.xy(xy, type, ...): "xlabl" is not a graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameter

Warning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameter
Warning in box(...): "xlabl" is not a graphical parameter
Warning in title(...): "xlabl" is not a graphical parameter

### 상호작용 효과(Interaction Effect)

res<-lm(S~X+E+M+E*M,data=data_5.1)

res

Call:
lm(formula = S ~ X + E + M + E * M, data = data_5.1)

Coefficients:
(Intercept)            X           E1           E2           M1        E1:M1  
    11203.4        497.0      -1730.7       -349.1       7047.4      -3066.0  
      E2:M1  
     1836.5  
summary(res)

Call:
lm(formula = S ~ X + E + M + E * M, data = data_5.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-928.13  -46.21   24.33   65.88  204.89 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11203.434     79.065 141.698  < 2e-16 ***
X             496.987      5.566  89.283  < 2e-16 ***
E1          -1730.748    105.334 -16.431  < 2e-16 ***
E2           -349.078     97.568  -3.578 0.000945 ***
M1           7047.412    102.589  68.695  < 2e-16 ***
E1:M1       -3066.035    149.330 -20.532  < 2e-16 ***
E2:M1        1836.488    131.167  14.001  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 173.8 on 39 degrees of freedom
Multiple R-squared:  0.9988,    Adjusted R-squared:  0.9986 
F-statistic:  5517 on 6 and 39 DF,  p-value: < 2.2e-16
### 표준화잔차 대 경력연수

plot(data_5.1$X,rstandard(res),pch=19,xlab="범주",ylab="잔차")

### Cook의 거리

plot(res,4) #33번째 데이터만 제외하고 다시 회귀모형을 생성예정

### 상호작용 효과 - 관측개체 33 제외 

data_use<-data_5.1[-33,]

res<-lm(S~X+E+M+E*M,data=data_use)

res

Call:
lm(formula = S ~ X + E + M + E * M, data = data_use)

Coefficients:
(Intercept)            X           E1           E2           M1        E1:M1  
    11199.7        498.4      -1741.3       -357.0       7040.6      -3051.8  
      E2:M1  
     1997.5  
summary(res)

Call:
lm(formula = S ~ X + E + M + E * M, data = data_use)

Residuals:
     Min       1Q   Median       3Q      Max 
-112.884  -43.636   -5.036   46.622  128.480 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11199.714     30.533 366.802  < 2e-16 ***
X             498.418      2.152 231.640  < 2e-16 ***
E1          -1741.336     40.683 -42.803  < 2e-16 ***
E2           -357.042     37.681  -9.475 1.49e-11 ***
M1           7040.580     39.619 177.707  < 2e-16 ***
E1:M1       -3051.763     57.674 -52.914  < 2e-16 ***
E2:M1        1997.531     51.785  38.574  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 67.12 on 38 degrees of freedom
Multiple R-squared:  0.9998,    Adjusted R-squared:  0.9998 
F-statistic: 3.543e+04 on 6 and 38 DF,  p-value: < 2.2e-16
### 표준화잔차 대 경력연수

plot(data_use$X,rstandard(res),pch=19,xlab="범주",ylab="잔차")

#### 표준화잔차 대 교육수준 - 관리 조합

EM<-paste0(data_use$E,data_use$M)

plot(EM,rstandard(res),pch=19,xlabl="범주",ylab="잔차")
Warning in plot.window(...): "xlabl" is not a graphical parameter
Warning in plot.xy(xy, type, ...): "xlabl" is not a graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameter

Warning in axis(side = side, at = at, labels = labels, ...): "xlabl" is not a
graphical parameter
Warning in box(...): "xlabl" is not a graphical parameter
Warning in title(...): "xlabl" is not a graphical parameter

#상호 호과가 들어간 이 모형이 더 괜찮은 모양이라고 판단이 된다 

### 기본 급료 추정 - 표 5.6

res<-lm(S~X+E+M+E*M,data=data_use)

df_new<-data.frame(X=rep(0,6),

                   E=rep(1:3,c(2,2,2)),

                   M=rep(c(0,1),3))

### 가변수 생성 - 분석용

df_new$E<-factor(as.character(df_new$E),levels = c("3","1","2")) 

df_new$M<-factor(as.character(df_new$M),levels = c("0","1")) 

cbind(df_new,predict = predict(res,df_new,interval = "confidence"))
  X E M predict.fit predict.lwr predict.upr
1 0 1 0    9458.378    9395.539    9521.216
2 0 1 1   13447.195   13382.933   13511.456
3 0 2 0   10842.672   10789.719   10895.624
4 0 2 1   19880.782   19814.090   19947.474
5 0 3 0   11199.714   11137.902   11261.525
6 0 3 1   18240.294   18182.503   18298.084
# 이를 보고 각 레벨에 따른 차이를 보고 얼마나 나는 지 분석이 가능해야 한다 

# ex) 고졸과 대학원졸의 관리자 직급의 급여의 차이는?(평균적으로)

5.4 회귀방정식의 체계: 두 집단의 비교

표 5.7 - 고용전 검사 프로그램 데이터

data_5.7<-read.table("All_Data/p140.txt",header=T,sep="\t")

head(data_5.7)
  TEST RACE JPERF
1 0.28    1  1.83
2 0.97    1  4.59
3 1.25    1  2.97
4 2.46    1  8.14
5 2.51    1  8.00
6 1.17    1  3.30
# 모형 1 - 통합모형 인종간 차이가 없을때

model_1<-lm(JPERF~TEST,data=data_5.7)

summary(model_1) 

Call:
lm(formula = JPERF ~ TEST, data = data_5.7)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3558 -0.8798 -0.1897  1.2735  2.3312 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.0350     0.8680   1.192 0.248617    
TEST          2.3605     0.5381   4.387 0.000356 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.591 on 18 degrees of freedom
Multiple R-squared:  0.5167,    Adjusted R-squared:  0.4899 
F-statistic: 19.25 on 1 and 18 DF,  p-value: 0.0003555
# 모형 3 

model_3<-lm(JPERF~TEST+RACE+TEST*RACE,data=data_5.7)

summary(model_3)

Call:
lm(formula = JPERF ~ TEST + RACE + TEST * RACE, data = data_5.7)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0734 -1.0594 -0.2548  1.2830  2.1980 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   2.0103     1.0501   1.914   0.0736 .
TEST          1.3134     0.6704   1.959   0.0677 .
RACE         -1.9132     1.5403  -1.242   0.2321  
TEST:RACE     1.9975     0.9544   2.093   0.0527 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.407 on 16 degrees of freedom
Multiple R-squared:  0.6643,    Adjusted R-squared:  0.6013 
F-statistic: 10.55 on 3 and 16 DF,  p-value: 0.0004511
# 인종적인 차이가 있는지 없는지를 확인해야 한다 어떤 모형을 사용할지 

### H_0:gamma=delta=0

anova(model_1,model_3) #model_3가 FM(완전모형)
Analysis of Variance Table

Model 1: JPERF ~ TEST
Model 2: JPERF ~ TEST + RACE + TEST * RACE
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     18 45.568                              
2     16 31.655  2    13.913 3.5161 0.05424 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#P-value(0.05424)<0.05이기에 H_0 이기에 모형1을 선택하는 것이 옳다고 판단(그러나 확신x)

#서로의 R-squared롤 보면 model_3가 더 좋음 / ANOVA는 참고용 절대적이지는 않다 

그림 5.7 - 표준화잔차 대 검사점수 : 모형 1

plot(data_5.7$TEST,rstandard(model_1),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.7 - 표준화잔차 대 검사점수 : 모형 1")

그림 5.8 - 표준화잔차 대 검사점수 : 모형 3

plot(data_5.7$TEST,rstandard(model_3),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.8 - 표준화잔차 대 검사점수 : 모형 3")

# 통계에서 나오는 결과는 결정의 보조수단이지 절대적이지 않아 

# 3번 모형을 선택한다고 결정한다고 진행 

summary(model_3) # Multiple R-squared:  0.6643

Call:
lm(formula = JPERF ~ TEST + RACE + TEST * RACE, data = data_5.7)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0734 -1.0594 -0.2548  1.2830  2.1980 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   2.0103     1.0501   1.914   0.0736 .
TEST          1.3134     0.6704   1.959   0.0677 .
RACE         -1.9132     1.5403  -1.242   0.2321  
TEST:RACE     1.9975     0.9544   2.093   0.0527 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.407 on 16 degrees of freedom
Multiple R-squared:  0.6643,    Adjusted R-squared:  0.6013 
F-statistic: 10.55 on 3 and 16 DF,  p-value: 0.0004511

그림 5.9 - 표준화잔차 대 검사점수 : 모형 1

plot(data_5.7$RACE,rstandard(model_1),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.9 - 표준화잔차 대 검사점수 : 모형 1")

# 분리된 회귀분석 결과

data_5.7_R1<-subset(data_5.7,RACE==1)

model_R1<-lm(JPERF~TEST,data=data_5.7_R1)

summary(model_R1) #소수민족

Call:
lm(formula = JPERF ~ TEST, data = data_5.7_R1)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0734 -0.6267 -0.2548  1.1624  1.5394 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.09712    1.03519   0.094 0.927564    
TEST         3.31095    0.62411   5.305 0.000724 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.292 on 8 degrees of freedom
Multiple R-squared:  0.7787,    Adjusted R-squared:  0.751 
F-statistic: 28.14 on 1 and 8 DF,  p-value: 0.0007239
data_5.7_R0<-subset(data_5.7,RACE==0)

model_R0<-lm(JPERF~TEST,data=data_5.7_R0)

summary(model_R0) #백인

Call:
lm(formula = JPERF ~ TEST, data = data_5.7_R0)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8599 -1.0663 -0.3061  1.0957  2.1980 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.0103     1.1291   1.780    0.113
TEST          1.3134     0.7208   1.822    0.106

Residual standard error: 1.512 on 8 degrees of freedom
Multiple R-squared:  0.2933,    Adjusted R-squared:  0.205 
F-statistic:  3.32 on 1 and 8 DF,  p-value: 0.1059
# 통합모형에서 나온 각각의 회귀식이 통합모형에서 나온것과 동일함 따라서 인종별로 나누어서

# 진행할 필요없이 통일모형을 사용해서 진행을 하면 된다(이는 데이터셋을 나눈경우와 동일함)

그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만

plot(data_5.7_R1$TEST,rstandard(model_R1),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만만")

그림 5.11 - 표준화잔차 대 검사점수 : 모형 1 . 백인만

plot(data_5.7_R0$TEST,rstandard(model_R0),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.11 - 표준화잔차 대 검사점수 : 모형 1. 백인만")

# 적절한 합격점수의 결정 - 소수민족

# 고용전 검사점수의 합격점에 대한 95% 신뢰구간

ym<-4

xm<-(ym-0.09712)/3.31095

s<-1.292

n<-10

t<-qt(1-0.05/2,8)

c(xm-(t*s/n)/3.31095, xm+(t*s/n)/3.31095) #신뢰구간
[1] 1.088795 1.268764

5.4.2 동일한 기울기와 다른 절편항을 가지는 모형

model_3<-lm(JPERF~TEST+RACE,data=data_5.7)

summary(model_3)

Call:
lm(formula = JPERF ~ TEST + RACE, data = data_5.7)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7872 -1.0370 -0.2095  0.9198  2.3645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.6120     0.8870   0.690 0.499578    
TEST          2.2988     0.5225   4.400 0.000391 ***
RACE          1.0276     0.6909   1.487 0.155246    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.54 on 17 degrees of freedom
Multiple R-squared:  0.5724,    Adjusted R-squared:  0.5221 
F-statistic: 11.38 on 2 and 17 DF,  p-value: 0.0007312
#Intercept = BETA_0 / TEST = BETA_1 / RACE = gamma

## H_0: gamma=0

anova(model_1,model_3) #p-value(0.1552)>0.05 이기에 H_0은 참이다//gamma=0
Analysis of Variance Table

Model 1: JPERF ~ TEST
Model 2: JPERF ~ TEST + RACE
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     18 45.568                           
2     17 40.322  1    5.2468 2.2121 0.1552
#기울기가 같고 절편이 다른 모형은 아니라고 데이터가 이야기하고 있다 

# 소수민족(RACE=1): (0.6120+1.0276)+2.2988*TEST = 1.6396+2.2988*TEST

# 백인(RACE=0): (0.6120)+2.2988*TEST

5.4.3 동일한 절편항과 다른 기울기를 가지는 모형

model_3<-lm(JPERF~TEST+I(TEST*RACE),data=data_5.7)

summary(model_3) #I()를 하면 교호작용하는 것만 보이게 하려고 없으면 RACE항이 자동추가됨

Call:
lm(formula = JPERF ~ TEST + I(TEST * RACE), data = data_5.7)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.41100 -0.88871 -0.03359  0.97720  2.44440 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)      1.1211     0.7804   1.437  0.16900   
TEST             1.8276     0.5356   3.412  0.00332 **
I(TEST * RACE)   0.9161     0.3972   2.306  0.03395 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.429 on 17 degrees of freedom
Multiple R-squared:  0.6319,    Adjusted R-squared:  0.5886 
F-statistic: 14.59 on 2 and 17 DF,  p-value: 0.0002045
## H_0: delta=0

anova(model_1,model_3) #p-value(0.03395)<0.05보다 작기에 delta항은 필요한 변수라는 사실 
Analysis of Variance Table

Model 1: JPERF ~ TEST
Model 2: JPERF ~ TEST + I(TEST * RACE)
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     18 45.568                              
2     17 34.708  1    10.861 5.3196 0.03395 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Intercept = BETA_0 / TEST = BETA_1 / I(TEST*RACE) = delta

6장 변수변환 : Transformation of Varibables

library(dplyr)

6.3 x-선 방사에 의한 박테리아 사망률

data_6.2<-read.table("All_Data/p168.txt",header=T,sep="\t")

head(data_6.2)
  t N_t
1 1 355
2 2 211
3 3 197
4 4 166
5 5 142
6 6 106

그림 6.5

plot(N_t~t,data=data_6.2,pch=19)

6.3.1 선형모형의 부적절성

res<-lm(N_t~t,data=data_6.2)

summary(res)

Call:
lm(formula = N_t ~ t, data = data_6.2)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.867 -23.599  -9.652  10.223 114.883 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   259.58      22.73  11.420 3.78e-08 ***
t             -19.46       2.50  -7.786 3.01e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 41.83 on 13 degrees of freedom
Multiple R-squared:  0.8234,    Adjusted R-squared:  0.8098 
F-statistic: 60.62 on 1 and 13 DF,  p-value: 3.006e-06

그림 6.6

plot(data_6.2$t,rstandard(res),pch=19) # 표준화 잔차의 플랏

# 잔차가 랜덤하게 잘 퍼져있어야 하는데 적절한 회귀모형이 아니라는 사실이 나옴

6.3.2 선형성을 위한 로그변환

그림 6.7

plot(log(N_t)~t,data=data_6.2,pch=19)

# 박테리아의 수에 로그를 취하니 선형성이 보인다 

res<-lm(log(N_t)~t,data=data_6.2)

summary(res)

Call:
lm(formula = log(N_t) ~ t, data = data_6.2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.18445 -0.06189  0.01253  0.05201  0.20021 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.973160   0.059778   99.92  < 2e-16 ***
t           -0.218425   0.006575  -33.22 5.86e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.11 on 13 degrees of freedom
Multiple R-squared:  0.9884,    Adjusted R-squared:  0.9875 
F-statistic:  1104 on 1 and 13 DF,  p-value: 5.86e-14
plot(data_6.2$t,rstandard(res),pch=19)

# 로그를 취해주니 잔차가 랜덤하게 이루어져 있음 

# n_0에 대한 추론

exp(5.973160)
[1] 392.7448
exp(coef(res)[1]) #로그를 취해주고 하는 부분이 잘 이해가 안됨 
(Intercept) 
   392.7449 
exp(coef(res)[1]-0.0588/2) #불편추정량 구하기 (참고용)
(Intercept) 
   381.3663 

6.4 분산안정화 변환

data_6.6<-read.table("All_Data/p174.txt",header=T,sep="\t")

data_6.6 #운항률과 사고 발생 건수에 대한 자료
   Y      N
1 11 0.0950
2  7 0.1920
3  7 0.0750
4 19 0.2078
5  9 0.1382
6  4 0.0540
7  3 0.1292
8  1 0.0503
9  3 0.0629
plot(Y~N,data=data_6.6,pch=19) #잔차의 분산이 계속 커지는 효과를 보임 

res_1<-lm(Y~N,data=data_6.6)

summary(res_1)

Call:
lm(formula = Y ~ N, data = data_6.6)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3351 -2.1281  0.1605  2.2670  5.6382 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -0.1402     3.1412  -0.045   0.9657  
N            64.9755    25.1959   2.579   0.0365 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.201 on 7 degrees of freedom
Multiple R-squared:  0.4872,    Adjusted R-squared:  0.4139 
F-statistic:  6.65 on 1 and 7 DF,  p-value: 0.03654
# 표준화잔차 대 N의 플롯 / 그림 6.11

plot(data_6.6$N,rstandard(res_1),pch=19) #등분산성에 만족하지 못하는 모형을 보인다 

# N에 대한 sqrt(Y)의 회귀

# N에 대한 Y의 회귀

res_2<-lm(sqrt(Y)~N,data=data_6.6)

summary(res_2)

Call:
lm(formula = sqrt(Y) ~ N, data = data_6.6)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9690 -0.7655  0.1906  0.5874  1.0211 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   1.1692     0.5783   2.022   0.0829 .
N            11.8564     4.6382   2.556   0.0378 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7733 on 7 degrees of freedom
Multiple R-squared:  0.4828,    Adjusted R-squared:  0.4089 
F-statistic: 6.535 on 1 and 7 DF,  p-value: 0.03776

그림 6.12

plot(data_6.6$N,rstandard(res_2),pch=19) #0을 중심으로 잔차가 보다 잘 퍼져있음이 보임 

6.5 이분산성의 검출

data_6.9<-read.table("All_Data/p176.txt",header=T,sep="\t")

data_6.9
      X   Y
1   294  30
2   247  32
3   267  37
4   358  44
5   423  47
6   311  49
7   450  56
8   534  62
9   438  68
10  697  78
11  688  80
12  630  84
13  709  88
14  627  97
15  615 100
16  999 109
17 1022 114
18 1015 117
19  700 106
20  850 128
21  980 130
22 1025 160
23 1021  97
24 1200 180
25 1250 112
26 1500 210
27 1650 135
# Y 대 X의 플로

plot(Y~X,data=data_6.9,pch=19,main="그림 6.13")

# X에 대한 Y의 회귀

res_1<-lm(Y~X,data=data_6.9)

summary(res_1)

Call:
lm(formula = Y ~ X, data = data_6.9)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.294  -9.298  -5.579  14.394  39.119 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.44806    9.56201   1.511    0.143    
X            0.10536    0.01133   9.303 1.35e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.73 on 25 degrees of freedom
Multiple R-squared:  0.7759,    Adjusted R-squared:  0.7669 
F-statistic: 86.54 on 1 and 25 DF,  p-value: 1.35e-09
# 표준화잔차 대 X의 플롯

plot(data_6.9$X, rstandard(res_1),pch=19,main = "그림 6.14")

6.6 이분산성의 제거

# 변환된 Y/X와 1/X를 적합한 회귀

data_6.9_1<-data.frame(Y=data_6.9$Y/data_6.9$X,

                       X=1/data_6.9$X)

res_2<-lm(Y~X,data=data_6.9_1)

summary(res_2) #지금은 변환하고 프라임 값들의 추정치가 나온것이기에 원래 회귀식으로 돌아가야함 

Call:
lm(formula = Y ~ X, data = data_6.9_1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.041477 -0.013852 -0.004998  0.024671  0.035427 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.120990   0.008999  13.445 6.04e-13 ***
X           3.803296   4.569745   0.832    0.413    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02266 on 25 degrees of freedom
Multiple R-squared:  0.02696,   Adjusted R-squared:  -0.01196 
F-statistic: 0.6927 on 1 and 25 DF,  p-value: 0.4131
# 본래 변환시 X를 나눴기에 X를 곱해준다 -> B_0,B_1의 추정값이 바뀐다 서로 

plot(data_6.9_1$X, rstandard(res_2),pch=19,

     xlab="1/X",ylab="잔차",main = "[그림 6.15]")

6.7 가중최소제곱법

wt<-1/data_6.9$X^2 #가중값

res_3<-lm(Y~X,data=data_6.9,weights = wt)

summary(res_3) #결과는 위의 6.6과 동일한 값이 나옴 

Call:
lm(formula = Y ~ X, data = data_6.9, weights = wt)

Weighted Residuals:
      Min        1Q    Median        3Q       Max 
-0.041477 -0.013852 -0.004998  0.024671  0.035427 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.803296   4.569745   0.832    0.413    
X           0.120990   0.008999  13.445 6.04e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02266 on 25 degrees of freedom
Multiple R-squared:  0.8785,    Adjusted R-squared:  0.8737 
F-statistic: 180.8 on 1 and 25 DF,  p-value: 6.044e-13

6.8 데이터에 대한 로그변환

data_6.9<-read.table("All_Data/p176.txt",header=T,sep="\t")

head(data_6.9)
    X  Y
1 294 30
2 247 32
3 267 37
4 358 44
5 423 47
6 311 49
# log(Y) 대 X의 산점도

plot((Y)~X,data=data_6.9,pch=19,main="식 6.9")

plot(log(Y)~X,data=data_6.9,pch=19,main="그림 6.16")

res_4<-lm(log(Y)~X,data=data_6.9)

summary(res_4)

Call:
lm(formula = log(Y) ~ X, data = data_6.9)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.59648 -0.16578  0.00244  0.17481  0.34964 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.5150232  0.1110670  31.648  < 2e-16 ***
X           0.0012041  0.0001316   9.153 1.85e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2524 on 25 degrees of freedom
Multiple R-squared:  0.7702,    Adjusted R-squared:  0.761 
F-statistic: 83.77 on 1 and 25 DF,  p-value: 1.855e-09
# X에 대한 log(Y)의 회귀로 부터 얻은 표준화잔차플롯

plot(data_6.9$X,rstandard(res_4),pch=19,

     xlab="X",ylab="잔차",main="그림 6.17")

# X와 X^2에 대한 log(Y)의 회귀

df<-data.frame(log_Y=log(data_6.9$Y),

               X=data_6.9$X,

               X2=data_6.9$X^2)

res_5<-lm(log_Y~X+X2,data=df) #그러나 이러한 과정은 필요업고 아래방식으로..

res_5<-lm(log(Y)~X+I(X^2),data=data_6.9)

summary(res_5)

Call:
lm(formula = log(Y) ~ X + I(X^2), data = data_6.9)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.30589 -0.11705 -0.02707  0.17593  0.30657 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.852e+00  1.566e-01  18.205 1.50e-15 ***
X            3.113e-03  3.989e-04   7.803 4.90e-08 ***
I(X^2)      -1.102e-06  2.238e-07  -4.925 5.03e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1817 on 24 degrees of freedom
Multiple R-squared:  0.8857,    Adjusted R-squared:  0.8762 
F-statistic: 92.98 on 2 and 24 DF,  p-value: 4.976e-12
# 표준화잔차 대 적합값 플롯 : X와 X^2에 대한 log(Y)의 회귀

plot(fitted(res_5),rstandard(res_5),pch=19,

     xlab="X",ylab="잔차",main="그림 6.18")

# 표준화잔차 대 X의 플롯 : X와 X^2에 대한 log(Y)의 회귀

plot(data_6.9$X,rstandard(res_5),pch=19,

     xlab="X",ylab="잔차",main="그림 6.19")

# 표준화잔차 대 X^2의 플롯 : X와 X^2에 대한 log(Y)의 회귀

plot(data_6.9$X^2,rstandard(res_5),pch=19,

     xlab="X",ylab="잔차",main="그림 6.20") #2차항을 추가함으로서 더 잘 피팅됨 

# 7장은 스킵 / 8장 스킵

9.2 통계적 추론에 미치는 효과

data_9.1<-read.table("All_Data/p236.txt",header=T,sep="\t")

head(data_9.1)
      ACHV      FAM     PEER   SCHOOL
1 -0.43148  0.60814  0.03509  0.16607
2  0.79969  0.79369  0.47924  0.53356
3 -0.92467 -0.82630 -0.61951 -0.78635
4 -2.19081 -1.25310 -1.21675 -1.04076
5 -2.84818  0.17399 -0.18517  0.14229
6 -0.66233  0.20246  0.12764  0.27311
res<-lm(ACHV~.,data_9.1)

summary(res)

Call:
lm(formula = ACHV ~ ., data = data_9.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2096 -1.3934 -0.2947  1.1415  4.5881 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06996    0.25064  -0.279    0.781
FAM          1.10126    1.41056   0.781    0.438
PEER         2.32206    1.48129   1.568    0.122
SCHOOL      -2.28100    2.22045  -1.027    0.308

Residual standard error: 2.07 on 66 degrees of freedom
Multiple R-squared:  0.2063,    Adjusted R-squared:  0.1702 
F-statistic: 5.717 on 3 and 66 DF,  p-value: 0.001535
#F-statistic: 5.717 on 3 and 66 DF, p-value: 0.001535

#이는 의미있는 beta가 존재한다라는 의미

#그러나 각 계수들의 회귀계수를 보니 모두 0이라는 결과가 나옴 이는 다중공선성이다

# Correlation panel

panel.cor<-function(x,y){

par(usr=c(0,1,0,1))

r<-round(cor(x,y),digits = 3)

text(0.5,0.5,r,cex=1.5)

}

pairs(data_9.1[-1],lower.panel = panel.cor)

그림 9.1

plot(fitted(res),rstandard(res),pch=19,

xlab="예측값",ylab="잔차",main="[그림 9.1]")

#이는 회귀모형에 동시에 들어가면 안된다라는 의미 (제거가 필요)

9.3 예측에 미치는 효과

data_9.5<-read.table("All_Data/p241.txt",header=T,sep="\t")

head(data_9.5)
  YEAR IMPORT DOPROD STOCK CONSUM
1   49   15.9  149.3   4.2  108.1
2   50   16.4  161.2   4.1  114.8
3   51   19.0  171.5   3.1  123.2
4   52   19.1  175.5   3.1  126.9
5   53   18.8  180.8   1.1  132.1
6   54   20.4  190.7   2.2  137.7
### 산점도

pairs(data_9.5[-1],lower.panel = panel.cor)

# 회귀분석(1) : 데이터 1949~1966

res<-lm(IMPORT~.,data_9.5[-1])

summary(res) #다중공선성이 존재한다고 예측을 일단함

Call:
lm(formula = IMPORT ~ ., data = data_9.5[-1])

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7208 -1.8354 -0.3479  1.2973  4.1008 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -19.7251     4.1253  -4.782 0.000293 ***
DOPROD        0.0322     0.1869   0.172 0.865650    
STOCK         0.4142     0.3223   1.285 0.219545    
CONSUM        0.2427     0.2854   0.851 0.409268    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.258 on 14 degrees of freedom
Multiple R-squared:  0.973, Adjusted R-squared:  0.9673 
F-statistic: 168.4 on 3 and 14 DF,  p-value: 3.212e-11

그림 9.3 : 표준화잔차의 인덱스플롯

plot(1:nrow(data_9.5),rstandard(res),

pch=19,type="b",

xla="번호",ylab="잔차",main="[그림 9.3]")

# 회귀분석(2) : 데이터 1949~1959

data_use<-subset(data_9.5,YEAR<=59)

res<-lm(IMPORT~.,data_use[-1])

summary(res)

Call:
lm(formula = IMPORT ~ ., data = data_use[-1])

Residuals:
     Min       1Q   Median       3Q      Max 
-0.52367 -0.38953  0.05424  0.22644  0.78313 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -10.12799    1.21216  -8.355  6.9e-05 ***
DOPROD       -0.05140    0.07028  -0.731 0.488344    
STOCK         0.58695    0.09462   6.203 0.000444 ***
CONSUM        0.28685    0.10221   2.807 0.026277 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4889 on 7 degrees of freedom
Multiple R-squared:  0.9919,    Adjusted R-squared:  0.9884 
F-statistic: 285.6 on 3 and 7 DF,  p-value: 1.112e-07

그림 9.4 : 잔차플롯

plot(1:nrow(data_use),rstandard(res),

pch=19,type="b",

xla="번호",ylab="잔차",main="[그림 9.4]")

9.4 다중공선성의 탐색

# 분산확대인자

library(olsrr)

# 교육기회 균등(EEO)

res_9.1<-lm(ACHV~.,data_9.1)

summary(res_9.1)

Call:
lm(formula = ACHV ~ ., data = data_9.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2096 -1.3934 -0.2947  1.1415  4.5881 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06996    0.25064  -0.279    0.781
FAM          1.10126    1.41056   0.781    0.438
PEER         2.32206    1.48129   1.568    0.122
SCHOOL      -2.28100    2.22045  -1.027    0.308

Residual standard error: 2.07 on 66 degrees of freedom
Multiple R-squared:  0.2063,    Adjusted R-squared:  0.1702 
F-statistic: 5.717 on 3 and 66 DF,  p-value: 0.001535
#car::vif(res_9.1)

olsrr::ols_vif_tol(res_9.1) # 10보다 크면 유의한 영향을 준다 제거필요?
  Variables  Tolerance      VIF
1       FAM 0.02660945 37.58064
2      PEER 0.03309981 30.21166
3    SCHOOL 0.01202567 83.15544

표 9.5 : 프랑스 경제 데이터

data_use<-subset(data_9.5,YEAR<=59)

res_9.5<-lm(IMPORT~.,data_use[-1])

#car::vif(res_9.1)

olsrr::ols_vif_tol(res_9.5) #작으면 새로운 변수를 만들거나 제거를 통해서 진행
  Variables   Tolerance        VIF
1    DOPROD 0.005376417 185.997470
2     STOCK 0.981441657   1.018909
3    CONSUM 0.005373166 186.110015
# 상태 지수

cnd.idx<-olsrr::ols_eigen_cindex(res_9.1)

round(cnd.idx,4)
  Eigenvalue Condition Index intercept    FAM   PEER SCHOOL
1     2.9547          1.0000    0.0005 0.0030 0.0037 0.0014
2     0.9974          1.7211    0.9756 0.0000 0.0000 0.0000
3     0.0400          8.5996    0.0004 0.3068 0.4428 0.0008
4     0.0079         19.2826    0.0235 0.6903 0.5535 0.9978
cnd.idx<-olsrr::ols_eigen_cindex(res_9.5)

round(cnd.idx,4)
  Eigenvalue Condition Index intercept DOPROD  STOCK CONSUM
1     3.8384          1.0000    0.0010 0.0000 0.0109 0.0000
2     0.1484          5.0863    0.0053 0.0001 0.9385 0.0001
3     0.0132         17.0732    0.7743 0.0015 0.0330 0.0011
4     0.0001        265.4613    0.2193 0.9984 0.0175 0.9989

10장 - 공선형 데이터의 처리

data_9.5<-read.table("All_Data/p241.txt",header=T,sep="\t")

head(data_9.5)
  YEAR IMPORT DOPROD STOCK CONSUM
1   49   15.9  149.3   4.2  108.1
2   50   16.4  161.2   4.1  114.8
3   51   19.0  171.5   3.1  123.2
4   52   19.1  175.5   3.1  126.9
5   53   18.8  180.8   1.1  132.1
6   54   20.4  190.7   2.2  137.7
# 회귀분석(2) : 데이터 1949~1959

data_use<-subset(data_9.5,YEAR<=59)

res<-lm(IMPORT~.,data_use[-1]) #1열 제거 

summary(res)

Call:
lm(formula = IMPORT ~ ., data = data_use[-1])

Residuals:
     Min       1Q   Median       3Q      Max 
-0.52367 -0.38953  0.05424  0.22644  0.78313 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -10.12799    1.21216  -8.355  6.9e-05 ***
DOPROD       -0.05140    0.07028  -0.731 0.488344    
STOCK         0.58695    0.09462   6.203 0.000444 ***
CONSUM        0.28685    0.10221   2.807 0.026277 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4889 on 7 degrees of freedom
Multiple R-squared:  0.9919,    Adjusted R-squared:  0.9884 
F-statistic: 285.6 on 3 and 7 DF,  p-value: 1.112e-07
head(data_use[-c(1:2)])
  DOPROD STOCK CONSUM
1  149.3   4.2  108.1
2  161.2   4.1  114.8
3  171.5   3.1  123.2
4  175.5   3.1  126.9
5  180.8   1.1  132.1
6  190.7   2.2  137.7
### 주성분(principle component)

pc<-prcomp(data_use[-c(1:2)],scale.=T)

pc$rotation
              PC1         PC2          PC3
DOPROD 0.70633041 -0.03568867 -0.706982083
STOCK  0.04350059  0.99902908 -0.006970795
CONSUM 0.70654444 -0.02583046  0.707197102
pc$x #변화된 새로운 변수가 저장된곳
          PC1         PC2         PC3
1  -2.1258872  0.63865815 -0.02072230
2  -1.6189273  0.55553922 -0.07111317
3  -1.1151675 -0.07297970 -0.02173008
4  -0.8942966 -0.08236998  0.01081318
5  -0.6442081 -1.30668523  0.07258248
6  -0.1903514 -0.65914745  0.02655252
7   0.3596219 -0.74367447  0.04278124
8   0.9718018  1.35405877  0.06286252
9   1.5593159  0.96404558  0.02357446
10  1.7669951  1.01521706 -0.04498818
11  1.9311034 -1.66266195 -0.08061267
### 주성분회귀 - 원래데이터를 주성분분석을 통해 새로운 데이터를 생성 이를 가지고 회귀 

df<-data.frame(IMPORT = scale(data_use$IMPORT),

               pc$x)

df
       IMPORT        PC1         PC2         PC3
1  -1.3185185 -2.1258872  0.63865815 -0.02072230
2  -1.2084753 -1.6189273  0.55553922 -0.07111317
3  -0.6362502 -1.1151675 -0.07297970 -0.02173008
4  -0.6142416 -0.8942966 -0.08236998  0.01081318
5  -0.6802675 -0.6442081 -1.30668523  0.07258248
6  -0.3281290 -0.1903514 -0.65914745  0.02655252
7   0.1780700  0.3596219 -0.74367447  0.04278124
8   1.0143989  0.9718018  1.35405877  0.06286252
9   1.3665374  1.5593159  0.96404558  0.02357446
10  1.2564942  1.7669951  1.01521706 -0.04498818
11  0.9703816  1.9311034 -1.66266195 -0.08061267
res<-lm(IMPORT~.,data=df)

summary(res) #유의한 1,2번째 계수들만 사용하고 3번째 것은 없이 해도 되겠다..

Call:
lm(formula = IMPORT ~ ., data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.11525 -0.08573  0.01194  0.04984  0.17236 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8.900e-16  3.244e-02   0.000 1.000000    
PC1         6.900e-01  2.406e-02  28.673 1.61e-08 ***
PC2         1.913e-01  3.406e-02   5.617 0.000801 ***
PC3         1.160e+00  6.559e-01   1.768 0.120376    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1076 on 7 degrees of freedom
Multiple R-squared:  0.9919,    Adjusted R-squared:  0.9884 
F-statistic: 285.6 on 3 and 7 DF,  p-value: 1.112e-07
#중간고사 이후의것이 주가나오겠지만 앞에것도 알아야함 연습문제에서 나올것 예상

11장

# 11.10 - 

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

head(data_3.3)
   Y X1 X2 X3 X4 X5 X6
1 43 51 30 39 61 92 45
2 63 64 51 54 63 73 47
3 71 70 68 69 76 86 48
4 61 63 45 47 54 84 35
5 81 78 56 66 71 83 47
6 43 55 49 44 54 49 34
# 분산확대 인자(VIF) : 10초과 > 심각한 공산성의 문제가 있음

res<-lm(Y~.,data=data_3.3)

summary(res)

Call:
lm(formula = Y ~ ., data = data_3.3)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9418  -4.3555   0.3158   5.5425  11.5990 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.78708   11.58926   0.931 0.361634    
X1           0.61319    0.16098   3.809 0.000903 ***
X2          -0.07305    0.13572  -0.538 0.595594    
X3           0.32033    0.16852   1.901 0.069925 .  
X4           0.08173    0.22148   0.369 0.715480    
X5           0.03838    0.14700   0.261 0.796334    
X6          -0.21706    0.17821  -1.218 0.235577    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared:  0.7326,    Adjusted R-squared:  0.6628 
F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05
olsrr::ols_vif_tol(res) #10을 넘는 것이 없기에 공산성의 문제는 없다 
  Variables Tolerance      VIF
1        X1 0.3749447 2.667060
2        X2 0.6246520 1.600891
3        X3 0.4403263 2.271043
4        X4 0.3248624 3.078226
5        X5 0.8142600 1.228109
6        X6 0.5124025 1.951591
# 수정결정계수

res_summ<-summary(res)

res_summ$adj.r.squared #0.662846
[1] 0.662846
# Mallow's Cp - 작을수록 좋음

# 축소모형

res_subset<-lm(Y~X1+X3,data=data_3.3)

olsrr::ols_mallows_cp(res_subset,res) 
[1] 1.114811
# AIC - 참고만

olsrr::ols_aic(res_subset,method = "SAS")
[1] 118.0024
# BIC

olsrr::ols_sbic(res_subset,res)
[1] 121.0938
### 전진적 선택법 - AIC

res_step<-step(res,direction = "forward")
Start:  AIC=123.36
Y ~ X1 + X2 + X3 + X4 + X5 + X6
### 후진적 제거법 -AIC

res_step<-step(res,direction = "backward")
Start:  AIC=123.36
Y ~ X1 + X2 + X3 + X4 + X5 + X6

       Df Sum of Sq    RSS    AIC
- X5    1      3.41 1152.4 121.45
- X4    1      6.80 1155.8 121.54
- X2    1     14.47 1163.5 121.74
- X6    1     74.11 1223.1 123.24
<none>              1149.0 123.36
- X3    1    180.50 1329.5 125.74
- X1    1    724.80 1873.8 136.04

Step:  AIC=121.45
Y ~ X1 + X2 + X3 + X4 + X6

       Df Sum of Sq    RSS    AIC
- X4    1     10.61 1163.0 119.73
- X2    1     14.16 1166.6 119.82
- X6    1     71.27 1223.7 121.25
<none>              1152.4 121.45
- X3    1    177.74 1330.1 123.75
- X1    1    724.70 1877.1 134.09

Step:  AIC=119.73
Y ~ X1 + X2 + X3 + X6

       Df Sum of Sq    RSS    AIC
- X2    1     16.10 1179.1 118.14
- X6    1     61.60 1224.6 119.28
<none>              1163.0 119.73
- X3    1    197.03 1360.0 122.42
- X1    1   1165.94 2328.9 138.56

Step:  AIC=118.14
Y ~ X1 + X3 + X6

       Df Sum of Sq    RSS    AIC
- X6    1     75.54 1254.7 118.00
<none>              1179.1 118.14
- X3    1    186.12 1365.2 120.54
- X1    1   1259.91 2439.0 137.94

Step:  AIC=118
Y ~ X1 + X3

       Df Sum of Sq    RSS    AIC
<none>              1254.7 118.00
- X3    1    114.73 1369.4 118.63
- X1    1   1370.91 2625.6 138.16
### 단계적 방법

res_step<-step(res) #최종적으로는 X1+X3이 선택됨 이거 사용하는 것이 좋을듯 
Start:  AIC=123.36
Y ~ X1 + X2 + X3 + X4 + X5 + X6

       Df Sum of Sq    RSS    AIC
- X5    1      3.41 1152.4 121.45
- X4    1      6.80 1155.8 121.54
- X2    1     14.47 1163.5 121.74
- X6    1     74.11 1223.1 123.24
<none>              1149.0 123.36
- X3    1    180.50 1329.5 125.74
- X1    1    724.80 1873.8 136.04

Step:  AIC=121.45
Y ~ X1 + X2 + X3 + X4 + X6

       Df Sum of Sq    RSS    AIC
- X4    1     10.61 1163.0 119.73
- X2    1     14.16 1166.6 119.82
- X6    1     71.27 1223.7 121.25
<none>              1152.4 121.45
- X3    1    177.74 1330.1 123.75
- X1    1    724.70 1877.1 134.09

Step:  AIC=119.73
Y ~ X1 + X2 + X3 + X6

       Df Sum of Sq    RSS    AIC
- X2    1     16.10 1179.1 118.14
- X6    1     61.60 1224.6 119.28
<none>              1163.0 119.73
- X3    1    197.03 1360.0 122.42
- X1    1   1165.94 2328.9 138.56

Step:  AIC=118.14
Y ~ X1 + X3 + X6

       Df Sum of Sq    RSS    AIC
- X6    1     75.54 1254.7 118.00
<none>              1179.1 118.14
- X3    1    186.12 1365.2 120.54
- X1    1   1259.91 2439.0 137.94

Step:  AIC=118
Y ~ X1 + X3

       Df Sum of Sq    RSS    AIC
<none>              1254.7 118.00
- X3    1    114.73 1369.4 118.63
- X1    1   1370.91 2625.6 138.16
summary(res_step) #X3가 의미 없다고 나와도 아님 검증된것임

Call:
lm(formula = Y ~ X1 + X3, data = data_3.3)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.5568  -5.7331   0.6701   6.5341  10.3610 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.8709     7.0612   1.398    0.174    
X1            0.6435     0.1185   5.432 9.57e-06 ***
X3            0.2112     0.1344   1.571    0.128    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.817 on 27 degrees of freedom
Multiple R-squared:  0.708, Adjusted R-squared:  0.6864 
F-statistic: 32.74 on 2 and 27 DF,  p-value: 6.058e-08

0613 - regression analysis

연습문제 11.5

data_use<-read.table("All_Data/p329.txt",header=T,sep="\t")

head(data_use)
     X1 X2    X3    X4 X5 X6 X7 X8 X9    Y
1 4.918  1 3.472 0.998  1  7  4 42  0 25.9
2 5.021  1 3.531 1.500  2  7  4 62  0 29.5
3 4.543  1 2.275 1.175  1  6  3 40  0 27.9
4 4.557  1 4.050 1.232  1  6  3 54  0 25.9
5 5.060  1 4.455 1.121  1  6  3 42  0 29.9
6 3.891  1 4.455 0.988  1  6  3 56  0 29.9
### 데이터 탐색 - 자료형

sapply(data_use,class)
       X1        X2        X3        X4        X5        X6        X7        X8 
"numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" 
       X9         Y 
"numeric" "numeric" 
### 데이터 탐색 - 기초 통계량

summary(data_use) #결측값의 여부를 확인
       X1              X2              X3              X4       
 Min.   :3.891   Min.   :1.000   Min.   :2.275   Min.   :0.975  
 1st Qu.:5.058   1st Qu.:1.000   1st Qu.:4.855   1st Qu.:1.161  
 Median :5.974   Median :1.000   Median :5.685   Median :1.432  
 Mean   :6.405   Mean   :1.167   Mean   :6.033   Mean   :1.384  
 3rd Qu.:7.873   3rd Qu.:1.500   3rd Qu.:7.158   3rd Qu.:1.577  
 Max.   :9.142   Max.   :1.500   Max.   :9.890   Max.   :1.831  
       X5              X6            X7              X8              X9      
 Min.   :0.000   Min.   :5.0   Min.   :2.000   Min.   : 3.00   Min.   :0.00  
 1st Qu.:1.000   1st Qu.:6.0   1st Qu.:3.000   1st Qu.:30.00   1st Qu.:0.00  
 Median :1.000   Median :6.0   Median :3.000   Median :40.00   Median :0.00  
 Mean   :1.312   Mean   :6.5   Mean   :3.167   Mean   :37.46   Mean   :0.25  
 3rd Qu.:2.000   3rd Qu.:7.0   3rd Qu.:3.250   3rd Qu.:48.50   3rd Qu.:0.25  
 Max.   :2.000   Max.   :8.0   Max.   :4.000   Max.   :62.00   Max.   :1.00  
       Y        
 Min.   :25.90  
 1st Qu.:29.90  
 Median :33.70  
 Mean   :34.63  
 3rd Qu.:38.15  
 Max.   :45.80  
### 데이터 탐색 - 히스토그램 & 상자그림

hist(data_use$Y,main="histogram of data_use$Y")

boxplot(data_use$Y)

### 데이터 탐색 - 산점도 행렬 & 상관계수

plot(data_use)

cor(data_use) 
           X1         X2         X3         X4          X5        X6        X7
X1  1.0000000  0.6512669  0.6892117  0.7342737  0.45855650 0.6406157 0.3671126
X2  0.6512669  1.0000000  0.4129558  0.7285916  0.22402204 0.5103104 0.4264014
X3  0.6892117  0.4129558  1.0000000  0.5715520  0.20466375 0.3921244 0.1516093
X4  0.7342737  0.7285916  0.5715520  1.0000000  0.35888351 0.6788606 0.5743353
X5  0.4585565  0.2240220  0.2046638  0.3588835  1.00000000 0.5893871 0.5412988
X6  0.6406157  0.5103104  0.3921244  0.6788606  0.58938707 1.0000000 0.8703883
X7  0.3671126  0.4264014  0.1516093  0.5743353  0.54129880 0.8703883 1.0000000
X8 -0.4371012 -0.1007485 -0.3527514 -0.1390869 -0.02016883 0.1242663 0.3135114
X9  0.1466825  0.2041241  0.3059946  0.1065612  0.10161846 0.2222222 0.0000000
Y   0.8739117  0.7097771  0.6476364  0.7077656  0.46146792 0.5284436 0.2815200
            X8        X9          Y
X1 -0.43710116 0.1466825  0.8739117
X2 -0.10074847 0.2041241  0.7097771
X3 -0.35275139 0.3059946  0.6476364
X4 -0.13908686 0.1065612  0.7077656
X5 -0.02016883 0.1016185  0.4614679
X6  0.12426629 0.2222222  0.5284436
X7  0.31351144 0.0000000  0.2815200
X8  1.00000000 0.2257796 -0.3974034
X9  0.22577960 1.0000000  0.2668783
Y  -0.39740338 0.2668783  1.0000000
pairs(data_use)

panel.cor<-function(x,y){

  par(usr=c(0,1,0,1))

  r<-round(cor(x,y),digits = 3)

  text(0.5,0.5,r,cex=1.5)

}

pairs(data_use,lower.panel = panel.cor)

### (a) 모든 변수들이 모형에 포함시킬 것인가?

cor(data_use)
           X1         X2         X3         X4          X5        X6        X7
X1  1.0000000  0.6512669  0.6892117  0.7342737  0.45855650 0.6406157 0.3671126
X2  0.6512669  1.0000000  0.4129558  0.7285916  0.22402204 0.5103104 0.4264014
X3  0.6892117  0.4129558  1.0000000  0.5715520  0.20466375 0.3921244 0.1516093
X4  0.7342737  0.7285916  0.5715520  1.0000000  0.35888351 0.6788606 0.5743353
X5  0.4585565  0.2240220  0.2046638  0.3588835  1.00000000 0.5893871 0.5412988
X6  0.6406157  0.5103104  0.3921244  0.6788606  0.58938707 1.0000000 0.8703883
X7  0.3671126  0.4264014  0.1516093  0.5743353  0.54129880 0.8703883 1.0000000
X8 -0.4371012 -0.1007485 -0.3527514 -0.1390869 -0.02016883 0.1242663 0.3135114
X9  0.1466825  0.2041241  0.3059946  0.1065612  0.10161846 0.2222222 0.0000000
Y   0.8739117  0.7097771  0.6476364  0.7077656  0.46146792 0.5284436 0.2815200
            X8        X9          Y
X1 -0.43710116 0.1466825  0.8739117
X2 -0.10074847 0.2041241  0.7097771
X3 -0.35275139 0.3059946  0.6476364
X4 -0.13908686 0.1065612  0.7077656
X5 -0.02016883 0.1016185  0.4614679
X6  0.12426629 0.2222222  0.5284436
X7  0.31351144 0.0000000  0.2815200
X8  1.00000000 0.2257796 -0.3974034
X9  0.22577960 1.0000000  0.2668783
Y  -0.39740338 0.2668783  1.0000000
model_1<-lm(data_use$Y~.,data_use)

summary(model_1) #전형적인 다중공선성이 보인다 

Call:
lm(formula = data_use$Y ~ ., data = data_use)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7729 -1.9801 -0.0868  1.6615  4.2618 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 15.31044    5.96093   2.568   0.0223 *
X1           1.95413    1.03833   1.882   0.0808 .
X2           6.84552    4.33529   1.579   0.1367  
X3           0.13761    0.49436   0.278   0.7848  
X4           2.78143    4.39482   0.633   0.5370  
X5           2.05076    1.38457   1.481   0.1607  
X6          -0.55590    2.39791  -0.232   0.8200  
X7          -1.24516    3.42293  -0.364   0.7215  
X8          -0.03800    0.06726  -0.565   0.5810  
X9           1.70446    1.95317   0.873   0.3976  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.973 on 14 degrees of freedom
Multiple R-squared:  0.8512,    Adjusted R-squared:  0.7555 
F-statistic: 8.898 on 9 and 14 DF,  p-value: 0.0002015
#coef가 X6 0.820이므로 확인 필요

olsrr::ols_vif_tol(model_1) #X6,X7이 10에 가깝기에 제거해야하는 대상중 최우선 
  Variables  Tolerance       VIF
1        X1 0.14241176  7.021892
2        X2 0.35267392  2.835480
3        X3 0.40735454  2.454864
4        X4 0.26066568  3.836332
5        X5 0.54842184  1.823414
6        X6 0.08539078 11.710866
7        X7 0.10286111  9.721847
8        X8 0.43083904  2.321052
9        X9 0.51482060  1.942424
#다중공선성이 보이기에 모든 변수를 모형에 포함 시킬수는 없다 

#10보다 크면 심각한 공선성이 존재 

### (b) 지방세(X1) 방의 수(X6), 건물의 나이(X8)가 판매가격(Y)을 

#       설명하는데 적절하다는 의견에 동의 하는가?

model_2<-lm(Y~X1+X6+X8,data=data_use)

summary(model_2)

Call:
lm(formula = Y ~ X1 + X6 + X8, data = data_use)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7486 -2.4082 -0.3594  2.1378  6.5353 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.796013   4.971105   2.976 0.007462 ** 
X1           3.489464   0.729368   4.784 0.000113 ***
X6          -0.415515   1.182262  -0.351 0.728921    
X8           0.004923   0.063597   0.077 0.939062    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.123 on 20 degrees of freedom
Multiple R-squared:  0.7655,    Adjusted R-squared:  0.7303 
F-statistic: 21.76 on 3 and 20 DF,  p-value: 1.653e-06
# VIF

olsrr::ols_vif_tol(model_2)
  Variables Tolerance      VIF
1        X1 0.3184367 3.140342
2        X6 0.3875669 2.580200
3        X8 0.5317389 1.880622
# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_2)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_2) 

#동의는 할 수 있으나(적절하지만) 최고의 모형이라는데는 동의 못함

### (c) 지방세 X1가 단독으로 판매가격 Y을 설명하는데 적절하다는 의견에 동의?

model_3<-lm(Y~X1,data=data_use)

summary(model_3) #adj-rsquared는 변수를 선택할때 사용 

Call:
lm(formula = Y ~ X1, data = data_use)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8445 -2.3340 -0.3841  1.9689  6.3005 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  13.3553     2.5955   5.146 3.71e-05 ***
X1            3.3215     0.3939   8.433 2.44e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.988 on 22 degrees of freedom
Multiple R-squared:  0.7637,    Adjusted R-squared:  0.753 
F-statistic: 71.11 on 1 and 22 DF,  p-value: 2.435e-08
# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_3)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_3) 

### 적절한 모형 제시

data_use_2<-data_use[-6]

model_4<-lm(Y~.,data_use_2)

summary(model_4)

Call:
lm(formula = Y ~ ., data = data_use_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7340 -1.8513 -0.0154  1.5472  4.2113 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 14.40540    4.36026   3.304  0.00482 **
X1           1.81642    0.82431   2.204  0.04360 * 
X2           7.13892    4.01353   1.779  0.09556 . 
X3           0.14721    0.47683   0.309  0.76177   
X4           2.73339    4.24921   0.643  0.52976   
X5           2.06520    1.33883   1.543  0.14377   
X7          -1.91236    1.79355  -1.066  0.30318   
X8          -0.03832    0.06509  -0.589  0.56481   
X9           1.48746    1.65931   0.896  0.38418   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.878 on 15 degrees of freedom
Multiple R-squared:  0.8506,    Adjusted R-squared:  0.771 
F-statistic: 10.68 on 8 and 15 DF,  p-value: 5.936e-05
# VIF

olsrr::ols_vif_tol(model_4)
  Variables Tolerance      VIF
1        X1 0.2117072 4.723504
2        X2 0.3855308 2.593827
3        X3 0.4102334 2.437637
4        X4 0.2612467 3.827800
5        X5 0.5495336 1.819725
6        X7 0.3510102 2.848920
7        X8 0.4310175 2.320091
8        X9 0.6683173 1.496295
# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_4)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_4) 

### 단계적 선택법 - AIC

model_5<-step(model_4)
Start:  AIC=57.45
Y ~ X1 + X2 + X3 + X4 + X5 + X7 + X8 + X9

       Df Sum of Sq    RSS    AIC
- X3    1     0.789 125.00 55.605
- X8    1     2.870 127.08 56.001
- X4    1     3.426 127.63 56.106
- X9    1     6.654 130.86 56.706
- X7    1     9.414 133.62 57.207
<none>              124.21 57.453
- X5    1    19.702 143.91 58.987
- X2    1    26.198 150.40 60.046
- X1    1    40.207 164.41 62.184

Step:  AIC=55.61
Y ~ X1 + X2 + X4 + X5 + X7 + X8 + X9

       Df Sum of Sq    RSS    AIC
- X8    1     3.369 128.36 54.243
- X4    1     4.816 129.81 54.513
- X9    1     9.581 134.58 55.378
- X7    1     9.655 134.65 55.391
<none>              125.00 55.605
- X5    1    18.957 143.95 56.994
- X2    1    25.474 150.47 58.057
- X1    1    53.245 178.24 62.122

Step:  AIC=54.24
Y ~ X1 + X2 + X4 + X5 + X7 + X9

       Df Sum of Sq    RSS    AIC
- X4    1     5.011 133.38 53.163
- X9    1     6.543 134.91 53.437
<none>              128.36 54.243
- X5    1    20.033 148.40 55.724
- X7    1    22.819 151.18 56.170
- X2    1    24.299 152.66 56.404
- X1    1    95.134 223.50 65.552

Step:  AIC=53.16
Y ~ X1 + X2 + X5 + X7 + X9

       Df Sum of Sq    RSS    AIC
- X9    1     6.223 139.60 52.257
<none>              133.38 53.163
- X5    1    17.801 151.18 54.169
- X7    1    17.873 151.25 54.181
- X2    1    39.335 172.71 57.365
- X1    1   157.972 291.35 69.915

Step:  AIC=52.26
Y ~ X1 + X2 + X5 + X7

       Df Sum of Sq    RSS    AIC
<none>              139.60 52.257
- X5    1    20.836 160.43 53.596
- X7    1    21.669 161.27 53.720
- X2    1    47.409 187.01 57.274
- X1    1   156.606 296.20 68.312
summary(model_5)

Call:
lm(formula = Y ~ X1 + X2 + X5 + X7, data = data_use_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5605 -2.0856  0.0238  1.8580  3.8981 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  13.6212     3.6725   3.709 0.001489 ** 
X1            2.4123     0.5225   4.617 0.000188 ***
X2            8.4589     3.3300   2.540 0.019970 *  
X5            2.0604     1.2235   1.684 0.108541    
X7           -2.2154     1.2901  -1.717 0.102176    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.711 on 19 degrees of freedom
Multiple R-squared:  0.8321,    Adjusted R-squared:  0.7968 
F-statistic: 23.54 on 4 and 19 DF,  p-value: 3.866e-07
# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_5)

layout(1)

# Cook의 거리 - 0.5보다 작으면 괜찮음 / 1보다 큰 값을 고려 / 전체적으로 봤을때 튀는 값이 있는 경우 

olsrr::ols_plot_cooksd_chart(model_5) 

### 17번 제거

data_use_3<-model_5$model[-17,]

model_6<-lm(Y~.,data_use_3)

summary(model_6)

Call:
lm(formula = Y ~ ., data = data_use_3)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4577 -1.6655  0.1575  1.7978  4.1865 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  13.2263     3.4659   3.816  0.00127 **
X1            1.7014     0.6247   2.723  0.01394 * 
X2           12.0705     3.6958   3.266  0.00429 **
X5            2.4602     1.1726   2.098  0.05028 . 
X7           -2.2310     1.2152  -1.836  0.08294 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.553 on 18 degrees of freedom
Multiple R-squared:  0.8418,    Adjusted R-squared:  0.8067 
F-statistic: 23.95 on 4 and 18 DF,  p-value: 5.316e-07
# VIF

olsrr::ols_vif_tol(model_6)
  Variables Tolerance      VIF
1        X1 0.3319061 3.012900
2        X2 0.3658998 2.732989
3        X5 0.5664655 1.765333
4        X7 0.6043771 1.654596
# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_6)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_6)

Source Code
---
title: "Regression Analysis"
author: "Hyunsoo Kim"
date: "2022-03-14"
categories: [Basic, Code, R]
page-layout: full
output:
  prettydoc::html_pretty:
    theme: architect
    highlight: github
editor_options: 
  chunk_output_type: console
mainfont: NanumGothic
---

# 0316 - 선형모형

> 예제를 통한 회귀분석

```{r}
getwd() #"C:/Users/Hyunsoo Kim/Documents/lecture/regression_analysis"
```

## 2장

### 표2.5

```{r}
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

dim(data_2.5) #14 2

head(data_2.5)

X<-data_2.5$Units

Y<-data_2.5$Minutes
```

### 표2.6

```{r}
df<-data.frame(

  #1:length(X),

  Y,

  X,

  Y-mean(Y),

  X-mean(X),

  (Y-mean(Y))^2,

  (X-mean(X))^2,

  (Y-mean(Y))*(X-mean(X))

)

df
```

### 공분산(Covariance) - 식2.2

```{r}
COV_XY<-sum((Y-mean(Y))*(X-mean(X))) / (length(X)-1) #136

### cov() 함수

cov(X,Y) #136

### 상관계수(correalationship)

### cor() 함수

cor(X,Y) #0.9936987 
```

## 0321 - 선형모형

```{r}
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

x<-data_2.5$Units

y<-data_2.5$Minutes
```

### 식 2.6

```{r}
cor_xy<- COV_XY / (sd(x)*sd(y))

cor_xy

### cor() 함수

cor(x,y)

cor(y,x)

data_2.5

cor(data_2.5)
```

### 그림 2.4

```{r}
class(X)

class(Y) #both numeric

plot(X,Y, pch=19,xlab="Units",ylab="Minutes") 
```

### 표 2.3

```{r}
data_2.3<-read.table("All_Data/p029a.txt",header=TRUE,sep="\t")

data_2.3

X<-data_2.3$X

Y<-data_2.3$Y
```

### 그림 2.2

```{r}
plot(X,Y)

cor(X,Y) # 0 완벽하게 2차함수의 형태도 0이 나옴(직선의 형태가 아닌것만)
```

### 표 2.4

```{r}
data_2.4<-read.table("All_Data/p029b.txt",header=TRUE,sep="\t")
```

### 그림 2.3

```{r}
plot(data_2.4$X1,data_2.4$Y1, pch=19); abline(3,0.5) #기울기 3 절편0.5인 선을 추가해라 

plot(data_2.4$X2,data_2.4$Y2, pch=19); abline(3,0.5)

plot(data_2.4$X3,data_2.4$Y3, pch=19); abline(3,0.5)

plot(data_2.4$X4,data_2.4$Y4, pch=19); abline(3,0.5)

m<-matrix(1:4,ncol=2,byrow=TRUE) #2행의 매트릭스 생성 

m

layout(m)

plot(Y1~X1, data=data_2.4,pch=19); abline(3,0.5)

#y~x  -> y=ax+b 이러한 형태를 가지는 모형식이라는 의미 

plot(Y2~X2, data=data_2.4,pch=19); abline(3,0.5)

plot(Y3~X3, data=data_2.4,pch=19); abline(3,0.5)

plot(Y4~X4, data=data_2.4,pch=19); abline(3,0.5)

layout(1) #변환을 다시 하지 않으면 설정한 매트릭스의 비율로 그래프가 그려짐 해제 필요 

# cor()

cor(data_2.4$X1,data_2.4$Y1) #0.8164205

cor(data_2.4$X2,data_2.4$Y2) #0.8162365

cor(data_2.4$X3,data_2.4$Y3) #0.8162867

cor(data_2.4$X4,data_2.4$Y4) #0.8165214

cor(data_2.4) #이렇게 한번에 할 수 있으나 가독성 떨어짐 
```

## 단순선형회귀모형 2.4

### 2.5 모수에 대한 추정

```{r}
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

x<-data_2.5$Units

y<-data_2.5$Minutes
```

### 식 2.14 & 2.15

```{r}
sum((y-mean(y))*(x-mean(x))) #1768

sum((x-mean(x))^2) #114

beta1_hat<-sum((y-mean(y))*(x-mean(x))) / sum((x-mean(x))^2)

beta1_hat #15.50877

beta0_hat <- mean(y) - (beta1_hat*mean(x))

beta0_hat #4.161654

### 최소제곱회귀 방정식

# Minutes = 4.161654 + 15.50877 * Units

plot(beta0_hat+beta1_hat*x,pch=19);

# 4개의 고장 난 부품을 수리하는데 걸리는 에측시간

4.161654 + 15.50877 * 4 #66.19673

units<-4

beta0_hat + beta1_hat*units

### 적합값(Fitted value)

y_hat<-beta0_hat + beta1_hat*(x)

### 최소 제곱 잔차(residual)

e<-y-y_hat

e #합이 0이라는 특징이 존재

sum(e) #1.278977e-13 0에 근사한 추지가 나옴
```

### 표2.7

```{r}
df_2.7<-data.frame(

  x=x,

  y=y,

  y_hat,

  e

)

df_2.7

### lm() 함수 (linear model)

# Minutes = beta0 + beta1 * Units + epsilon

# 모형식 : y~x

lm(y~x)

res_lm<-lm(Minutes~Units,data=data_2.5)

res_lm

# 리스트의 이름 

names(res_lm)

# 회귀계수

res_lm$coefficients

coef(res_lm)

# 적합값

res_lm$fitted.values

fitted(res_lm)

# 최소제곱잔차

res_lm$residuals

resid(res_lm)

residuals(res_lm)
```

### 그림 2.5 - 산점도 + 회귀선

```{r}
plot(beta0_hat+beta1_hat*x,pch=19);

#abline(beta0_hat,beta1_hat)

abline(res_lm)
```

## 0323

```{r}
data_2.5<-read.table("All_Data/p031.txt",header=TRUE,sep="\t")

res_lm <- lm(Minutes ~ Units, data = data_2.5)

res_lm

res_lm_summ<-summary(res_lm)

res_lm_summ #Pr(>|t|) - p-value

# unit은 시간에 영향을준다 약15.5분 만큼씩 

# coefficient에서 p-value에 대해서 알 수 있음 

# beta_0는 0이라고 보면되느냐? p-value가 0.05보다 크기에 
```

## 0328

### 2.7 신뢰구간

```{r}
confint(res_lm) # beta_0,1의 95% 신뢰구간을 뽑아줌 

?confint #level = 1-alpha

confint(res_lm, level=0.90) # 90%의 신뢰구간
```

### 2.8 예측

```{r}
# 4개의 고장난 부품을 수리하는 데에 걸리는 시간 예측

x<-4

4.161654 + 15.508772 *4

res_lm$coefficients[1]+(res_lm$coefficients[2]*x)

### predict()

df<-data.frame(Units=4) 

predict(res_lm,newdata=df) # res_lm을 만들때 사용한 데이터형식으로 만들어주어야함

res_lm_pred<-predict(res_lm,newdata=df,se.fit=TRUE)

### 예측값

res_lm_pred$fit

### 표준오차

res_lm_pred$se.fit # 평균반응에 대한 표준오차 

### 예측한계

df<-data.frame(Units=4) #예제서는 4대기준

df<-data.frame(Units=1:10) #다른것도 보고 싶은경우 

res_lm_pred_int_p<-predict(res_lm,newdata=df,interval="prediction")

### 신뢰한계

df<-data.frame(Units=1:10) #다른것도 보고 싶은경우 

res_lm_pred_int_c<-predict(res_lm,newdata=df,interval="confidence") #둘의 차이를 보면 예측한계의 범위가 더큼 

### 예측한계 & 신뢰한계

# 신뢰한계는 평균에서 멀어지만 오차의범위가 커지고 평균에 다가갈수록 오차가 줄어듬

plot(Minutes~Units,data=data_2.5,pch=19)

abline(res_lm,col="red",lwd=2)

lines(1:10,res_lm_pred_int_p[,"lwr"],col="darkgreen")

lines(1:10,res_lm_pred_int_p[,"upr"],col="darkgreen")

lines(1:10,res_lm_pred_int_c[,"lwr"],col="blue")

lines(1:10,res_lm_pred_int_c[,"upr"],col="blue")
```

### 2.9 적합성의 측정

```{r}
res_lm_summ<-summary(res_lm)

res_lm_summ #Multiple R-squared:0.9874 -> 반응변수의 전체변이중 98.94%가 예측변수에 의해 설명된다

# 만약 R-squared가 1이면 완벽한 선형의 관계 100%라는 것을 의미한다.

# R-squared는 변수가 들어갈수록 커지기에 adjust R-squared를 사용 추후 설명 
```

### 2.10 원점을 통과하는 회귀선

```{r}
# Minutes = beta1 + Units + epsilon

res_lm_no<-lm(Minutes~Units-1,data=data_2.5)

res_lm_no

summary(res_lm_no)

coef(summary(res_lm_no)) #rsquared=0.9975
```

### 2.11

```{r}
y<-rnorm(30)

t.test(y,mu=0)

summary(lm(y~1))
```

## 3장 다중선형회귀

### 3.3 사례 감독자 직무수행능력 데이터

```{r}
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

dim(data_3.3)

class(data_3.3)

sapply(data_3.3,class) #all numeric

summary(data_3.3) #모든변수가 numeric이면 분위수도 보여준다 

### 산점도 행렬

plot(data_3.3)
```

### 3.4 모수 추정

```{r}
lm(Y~X1+X2+X3+X4+X5+X6,data=data_3.3)

lm(Y~.,data=data_3.3) # X1+X2+X3+X4+X5+X6쓰는 것이 아니라 .을 써서 모든 변수를 다써줌 
```

### 3.5 회귀계수에 대한 해석

```{r}
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

lm(Y~X1+X2,data=data_3.3)

# (Intercept)         X1           X2  

#  15.32762      0.78034     -0.05016

# 1) Y에서 X1 효과 제거

m1<-lm(Y~X1,data=data_3.3) # y prime

m1$residuals # x1이 설명하지 못한값 / x1의 효과를 제거한 값

# 2) X2에서 X1 효과 제거

m2<-lm(X2~X1,data=data_3.3) # x2 prime 

m2$residuals 

# 3) X1의 효과가 제거된 Y와 X2의 적합 - 원점을 지나는 회귀선

lm(m1$residuals~m2$residuals-1) # 원점을 지나면 -1를 하고 진행 // -3.25e-17

# 다른 효과 없이(다른값이 고정) Y에 영향을 주는 순수한 X2의 값

# m2$residuals  : -0.05016  ==  X2 : -0.05016  

### 단위길이 척도화 - 잘사용하지않음

fn_scaling_len<-function(x){

  x0<-x-mean(x)

  x0/sqrt(sum(x0^2))

}

data_3.3_len<-sapply(data_3.3, fn_scaling_len)

data_3.3_len<-data.frame(data_3.3_len)

summary(data_3.3_len)

lm(Y~.,data=data_3.3_len)

### 표준화

# scale()

data_3.3_std<-scale(data_3.3)

#summary(data_3.3_std)

#sapply(data_3.3_std, sd, na.rm=T)

#class(data_3.3_std) #"matrix"

data_3.3_std<-data.frame(data_3.3_std)

#class(data_3.3_std) #"data.frame"

#sapply(data_3.3_std, sd, na.rm=T)

lm(Y~.,data=data_3.3_std) # beta게수 구하기 

```

### 3.7 최소제곱추정량의 성질

```{r}
res_lm<-lm(Y~.,data=data_3.3)

res_lm

summary(res_lm)

m1<-summary(lm(Y~X1+X2+X3+X4+X5+X6,data=data_3.3)) #Adjusted R-squared:  0.6628 

m2<-summary(lm(Y~X1+X2+X3+X4+X5,data=data_3.3)) #Adjusted R-squared:  0.6561 

# X6가 들어가는 것이 더 좋은 모델 

m1$adj.r.squared

m2$adj.r.squared #summary에서 보다 더 정확하게 수치가 나옴 
```

### 3.9 개별 회귀계수들에 대한 추론

```{r}
res_lm_summ<-summary(res_lm)

res_lm_summ #p-value의 존재는 무언가를 검정했다라는 반증

# p-value<0.05 H_1 귀무가설 채택 

# p-value>0.05 H_0 영가설 채택 // X1을 제외하고는 영가설 유의한 의미가 없음(Y에영향주는)

# 모두다 0이라는 가설을 가지고 분모 분자의 오차가 카이제곱을 따르고 거기서 나온 통계량

# F-분포 자유도는 분자 분모 두개를 가짐 //모아서 계산을 하기에 각각 계산하는것과 결과다름 

# 영가설-모든 회귀계수가 0이다.

# 대립가설-적어도 하나는 0이 아니다. p-value: 1.24e-05 <0.05 대립가설 채택 

# p-value가 0.05보다 작으면 대립가설 채택!!!!!! 기억해 

# 회귀계수에 대한 신뢰구간 - 95% 신뢰한계

confint(res_lm) #-13.18712881 ~ 34.7612816

#X1  0.28016866  0.9462066  사이에 0이 들어가있으면 영향을 준다라느걸 의미

#X2 -0.35381806  0.2077178  p-value없이도 알 수 있음 

#X5가 가장 영향이 적음 p-value가 가장 크기에(영향 효과의 크기를 비교할때)

#p-value가 작을 수록 영향을 많이 준다 beta값을 보는 것이 아닌 p-value를 보는 것 중요

#가장 의미있는 변수?->p-value가장 작은거 // 대립가설채택 Y에 영향을 가장

```

### 3.10 선형모형에서의 가설검정

### 3.10.1 모든 회귀예수들이 0인가에 대한 검정

```{r}
# H_0: beta_1:beta_6=0

model_reduce<-lm(Y~1,data=data_3.3)

model_full<-lm(Y~.,data=data_3.3)

anova(model_reduce,model_full)

#대립가설 = 완전모형이 적절하다 / 1.24e-05 *** < 0.05 

#의미 있는 예측 변수가 한개 이상 존재한다 

summary(model_full) #summary에서 beta_1~beta_6까지 모두가 0이라는 가설로 진행을 이미함

#F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05

#예상// 가장의미있는변수? -> X1 이유?-> p-value 0.000903 로 가장 작기에 영향많이 줄것으로 예측 
```

### 3.10.2 회귀계수들의 부분집합이 0인가에 대한 검정

```{r}
model_reduce<-lm(Y~X1+X3,data=data_3.3)

model_full<-lm(Y~.,data=data_3.3)

anova(model_reduce,model_full) #0.7158 > 0.05

#영가설은 H_0: b_2=b_4=b_5=b_6=0 이라는 사실을 알 수 있다 

#b_1&b_3는 반응변수에 유의한 반응을 준다라는 것도 연계하여 알 수 있다 
```

### 3.10.3 회귀계수들의 통일성에 대한 검정

```{r}
#해당 조건이 주어지고 만족할 때 beta_1=beta_3은 맞는가?

model_reduce<-lm(Y~I(X1-X3),data=data_3.3) #I를 씌우면 새로운 변수를 만든것과 동일

# X1-X3를 한 그자체를 분석하라는 의미//본래는 X1-X3 해서 새로운 변수를 만들어서 해야함 

model_full<-lm(Y~X1+X3,data=data_3.3)

anova(model_reduce,model_full) 

#install.packages("car")

library(car)

model_full<-lm(Y~X1+X3,data=data_3.3)

car::linearHypothesis(model_full,c("X1=X3"))
```

### 3.10.4 제약조건하에서 회귀계수에 대한 추정과 검정

```{r}
# H_0: beta_1+beta_3=1 | beta_2=beta_3:beta_6=0

model_full<-lm(Y~X1+X3,data=data_3.3)

car::linearHypothesis(model_full,c("X1 + X3 = 1"))

# x1의 효과가 증가하면 x3의 효과는 감소한다 상대적인 관계 (반대로도 가능)
```

### 3.11 예측

```{r}
model_full<-lm(Y~.,data_3.3)

# 예측값 - 적합값

model_full$fitted.values

# 예측한계(Prediction Limits)

predict(model_full,newdata = data_3.3,interval = "prediction")

# 신뢰한계(Confidence limits)

predict(model_full,newdata = data_3.3,interval = "confidence")
```

## 부록 : 행렬을 이용한 회귀계수 추정

```{r}
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

Y<-data_3.3$Y

X<-data_3.3[,-1]

X<-cbind(1,X)

X<-as.matrix(X)

#beta_hat<-solve(t(X) %*% X) %*% t(X) %*% Y # %*%행렬 계산 

P<-solve(t(X) %*% X) %*% t(X)

beta_hat<- P %*% Y

lm(Y~.,data_3.3)
```

## 4장 회귀진단: 모형위반의 검출

### 4.3 다양한 유형의 잔차들

```{r}
# 표준화 잔차

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

res_lm<-lm(Y~.,data=data_3.3)

class(res_lm)

mode(res_lm)

names(res_lm)

res_lm$fitted.values

str(res_lm)

# 잔차 

res_lm$residuals

resid(res_lm) #실제값에서 예측된 값을 뺸값

### 내적 표준화잔차

rstandard(res_lm)

### 외적 표준화잔차

MASS::studres(res_lm)

redsid_df<-data.frame(

  Y=data_3.3$Y,

  Y_hat=res_lm$fitted.values,

  resid=resid(res_lm),

  rstandard=rstandard(res_lm),

  studres=MASS::studres(res_lm)

)

redsid_df
```

### 4.5 모형을 적합깅 이전의 그래프

### 4.5.1 일차원 그래프

```{r}
a<-rnorm(100,70,10) #연속형 데이터

# 히스토그램 

hist(a)

hist(a,breaks=5) #범위를 조절 막대의 5번 자름 

# 줄기 잎 그림 

stem(a)

stem(round(a)) #줄기잎을 그릴때는 반올림을 하고 항상 진행 

stem(round(a),scale=2) #scale을 2배로 늘려라 5기준으로 반으로 잘라서 

# 모든데이터를 볼 수있는 장점 데이터가 많으면 구림 

# 점플롯

idx<-rep(1,length(a)) #a의 갯수에 맞춰서 1를 반복 

plot(idx,a)

plot(jitter(idx),a,xlim=c(0.5,1.5))

# 상자그림

boxplot(a) #사분위수에 대해서 알 수 있음 

# 상자그림 + 점플롯

boxplot(a)

points(jitter(idx),a)
```

### 4.5.2 이차원 그래프

```{r}
data_4.1<-read.table("All_Data/p103.txt",header=T,sep="\t")

data_4.1

class(data_4.1) #data.frame

# 산점도 행렬

plot(data_4.1)

cor(data_4.1) #상관계수

pairs(data_4.1)

# Correlation panel

panel.cor<-function(x,y){

  par(usr=c(0,1,0,1))

  r<-round(cor(x,y),digits = 3)

  text(0.5,0.5,r,cex=1.5)

}

pairs(data_4.1,lower.panel = panel.cor)

# 회전도표, 동적 그래프(3차원)

#install.packages("rgl")

library(rgl)

plot3d(x=data_4.1$X1,y=data_4.1$X2,z=data_4.1$Y) 
```

### 4.7 선형성과 정규성 가정에 대한 검토

```{r}
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

res_lm<-lm(Y~X1+X3,data=data_3.3)

layout(matrix(1:4,nrow=2,byrow=T))

plot(res_lm) #회귀진단 플랏이 나옴 

layout(1)

# 1. 표준화잔차의 정규확률분포

plot(res_lm,2)

# 2. 표준화잔차 대 각 예측변수들의 산점도

plot(res_lm,3) #이런경우에는 데이터를 추가한다 좌측하단이 없어서 

# 3. 표준화잔차 대 적합값의 플롯

# 4. 표준화잔차의 인덱스 플롯
```

## 0502 - 선형모형

```{r}
library(dplyr)
```

### 4.7 선형성과 정규성 가정에 대한 검토

```{r}
data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

res_lm<-lm(Y~X1+X3,data=data_3.3) #두개의 변수만 의미있다고 가정하고 진행 

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(res_lm) #회귀진단 플랏이 나옴 / 총6개임 

layout(1) #다시 한개의 플랏만 그리도록 

# 1. 표준화잔차의 정규확률분포

plot(res_lm,2) #QQ-plot y=x 기울기의 직선위에 점들이 있어야 한다 눈대중으로 

# 2. 표준화잔차 대 각 예측변수들의 산점도

plot(res_lm,3) #이런경우에는 데이터를 추가한다 좌측하단이 없어서 

#랜덤하게 데이터가 흩어져 있어야 한다

# 내적 표준화잔차

plot(data_3.3$X1,rstandard(res_lm))

plot(data_3.3$X3,rstandard(res_lm)) #각각의 잔차들이 랜덤하게 잘 퍼져야함 

# 3. 표준화잔차 대 적합값의 플롯

plot(res_lm,1) #잔차와 적합값은 상관성이 없어야하며 랜덤하게 퍼져야함 

# 4. 표준화잔차의 인덱스 플롯 

plot(res_lm,5)

data_2.4<-read.table("All_Data/p029b.txt",header=T,sep="\t")

m<-matrix(1:4,ncol=2,byrow=TRUE)  

layout(m)

plot(Y1~X1, data=data_2.4,pch=19); abline(3,0.5)

plot(Y2~X2, data=data_2.4,pch=19); abline(3,0.5)

plot(Y3~X3, data=data_2.4,pch=19); abline(3,0.5)

plot(Y4~X4, data=data_2.4,pch=19); abline(3,0.5)

layout(1)
```

### Figure4.4

```{r}
res_lm<-lm(Y1~X1,data=data_2.4) 

#1번플랏은 적당히 잘퍼짐,2번플랏은 어느정도 선형성 있음(데이터적어서그런거임)

res_lm<-lm(Y2~X2,data=data_2.4)

res_lm<-lm(Y3~X3,data=data_2.4)

res_lm<-lm(Y4~X4,data=data_2.4)

res_lm

layout(matrix(1:4,nrow=2,byrow=T))

plot(res_lm) 

layout(1)
```

### 4.8 지레점, 영향력, 특이값

```{r}
# 사례: 뉴욕 강 데이터

# agr-농업, forest-숲, rsdntial-주거, comlndl-산업, nitrogen-질소

data_1.9<-read.table("All_Data/p010.txt",header=T,sep="\t")

head(data_1.9)

plot(data_1.9[-1],pch=19) #river라는 첫번째 컬럼을 제외하고 진행 

res_1<-lm(Nitrogen~.,data=data_1.9[-1]) #모든데이터 사용

res_2<-lm(Nitrogen~.,data=data_1.9[-4,-1]) #4번쨰 데이터 제외

res_3<-lm(Nitrogen~.,data=data_1.9[-5,-1]) #5번쨰 데이터 제외 

#회귀계수

data.frame(all=coef(res_1),

           rm4=coef(res_2),

           rm5=coef(res_3))

#p-value

data.frame(all=coef(summary(res_1))[,4],

           rm4=coef(summary(res_2))[,4],

           rm5=coef(summary(res_3))[,4])

#4,5번째 데이터를 각각 뺴고 진행을 해보니 영향을 끼치는 값임을 알수있고

#5번째는 주거 관련해서는 부호를 바꿀 정도로 강력하다 

# 단순선형회구모형

res<-lm(Nitrogen~ComIndl,data=data_1.9)
```

### 그림 \[4.5\]

```{r}
plot(Nitrogen~ComIndl,data=data_1.9,pch=19)

abline(res)

#leverage values 지레값

p_ii<-hatvalues(res)

hiegh_leverage<-ifelse(p_ii>2*2/30,data_1.9$River,"")

hiegh_leverage #높은 지레값을 가지고 있는 강의 이름을 표시(이는 보기에 편하기 위해서함)

text(data_1.9$ComIndl,data_1.9$Nitrogen-0.1,hiegh_leverage)
```

### 그림 \[4.6\] - 잔차의 인덱스플롯

```{r}
plot(rstandard(res),pch=19) #2또는 3시그마를 넘으면 특이값이라 함 
```

### 그림 \[4.6\] - 지레값의 인덱스플롯

```{r}
plot(p_ii,pch=19) #평균의 2배를 기준으로 비교함 

abline(h=2*2/30,col="red") #이보다 높은것이 지레값이 높은것 높은 영향력을 가진것 

# 단순선형회귀모형

res<-lm(Nitrogen~ComIndl,data=data_1.9)

plot(Nitrogen~ComIndl,data=data_1.9,pch=19)

abline(res)

res<-lm(Nitrogen~ComIndl,data=data_1.9[-4,-1])

plot(Nitrogen~ComIndl,data=data_1.9[-4,-1],pch=19)

abline(res) #4번째 제외

res<-lm(Nitrogen~ComIndl,data=data_1.9[-5,-1])

plot(Nitrogen~ComIndl,data=data_1.9[-5,-1],pch=19)

abline(res) #5번쨰 제외

#이처럼 4,5번을 빼고 진행을 하면 조금더 잘 나타냄

res<-lm(Nitrogen~ComIndl,data=data_1.9[c(-4,-5),-1])

plot(Nitrogen~ComIndl,data=data_1.9[-5,-1],pch=19)

abline(res) #4,5번째 제외 
```

### 4.9 영향력의 측도

### 4.9.1 Cook의 거리

```{r}
res<-lm(Nitrogen~ComIndl,data=data_1.9)

plot(res,4)

#install.packages("olsrr")

library(olsrr)

ols_plot_cooksd_chart(res)
```

### 4.9.2 Welsch & Kuh의 측도

```{r}
olsrr::ols_plot_dffits(res)
```

### 4.9.3 Hadi의 영향력 측도

```{r}
olsrr::ols_plot_hadi(res)

# Residual & Leverage & Cook's distance

plot(res,5) #영향력 관측치를 보기 위한 플랏 

```

### 4.10 잠재성 - 잔차플롯

```{r}
olsrr::ols_plot_resid_pot(res) #2.0에 있는 값외에도 x축 0.2이상의 것들도 특이값으로 

#지레값이 커도 영향력이 없는 애들은 신경 안써도 되나 영향력이 큰애들을 탐색해 보아야 함
```

### 4.11 특이값에 대한 처리

```{r}
#특이값이 큰경우 그것이 TRUE데이터면 오류가 있는 경우 수정을 하거나 가중치를 변화를

#하거나 데이터를 수정을 시켜주거나 다시 실험을 하는 등 여러가지 방법을 사용하여.. 

## 첨가변수플롯(added-variable plot), 편회귀플롯(partial regression plot)

res<-lm(Nitrogen~ComIndl,data=data_1.9)

car::avPlots(res,pch=19)

res<-lm(Nitrogen~.,data=data_1.9[-1])

car::avPlots(res,pch=19)

#beta별로 각각의 어떤 변수가 영향력을 많이 주는지 알게하는 함수 
```

### 사례:스코틀랜드 언덕 경주 데이터

```{r}
data_4.5<-read.table("All_Data/p120.txt",header=T,sep="\t")

dim(data_4.5)

names(data_4.5)

head(data_4.5)

# 회전도표 

library(rgl)

with(data_4.5,plot3d(x=Distance,y=Climb,z=Time))

res<-lm(Time~Distance+Climb,data=data_4.5)

summary(res) #beta_0가 마이너스여도 신경안씀 관심있는 것은 회귀계수임 

#Time=-539.4829+373.0727*Distance+0.6629*Climb 
```

```{r}
### 첨가변수플롯(added-variable plot), 편회귀플롯(partial regression plot)

car::avPlots(res,pch=19)

### 성분잔차플롯(component plus residual plots), 편자차플롯(partial residual plot)

car::crPlots(res,id=T,pch=19) #첨가변수보다 성분잔차를 더 많이 사용 / 비선형여부를 확인

# 점선에 비해서 분홍선이 크게 떨어져 있지 않아 선형적인 추세를 가지고 있다고 추정이 가능 

### 잠재성-잔차플롯

olsrr::ols_plot_resid_pot(res)

### Hadi의 영향력 측도

olsrr::ols_plot_hadi(res)

### Cook의 거리

olsrr::ols_plot_cooksd_chart(res) #전체적은 플롯을 보면서 이상치에 대한 확인을 함 

# 이러한 값들의 제외 여부는 연구자가 선택해서 진행을 함 

# outlier에 대한 처리를 어떻게 했다고 말을 해야함 

# 보고서는 옆에 사람이 보고 쉽게 따라할 수 있을 정도로 

# 어떤 속성으로 어떻게 진행을 했다라는 것을 중시 코드 보단 결과를 보여줘라 
```

### 4.14 로버스트 회귀 - outlier에 크게 영향을 받지 않음

```{r}
library(MASS)

res<-lm(Time~Distance+Climb,data=data_4.5)

res

res_rlm<-MASS::rlm(Time~Distance+Climb,data=data_4.5)

res_rlm

summary(res)

summary(res_rlm)

#install.packages("robustbase")

library(robustbase)

res_lmrob<-lmrob(Time~Distance+Climb,data=data_4.5)

res_lmrob

summary(res_lmrob)

#standard error가 res > res_rlm > res_lmrob 순으로 되어 있음 

# 4장 연습문제 해보기
```

## 0509 - regression analysis

```{r}
library(dplyr)
```

## 5장 질적 예측 변수

```{r}
#install.packages("fastDummies")

library(fastDummies)
```

## 5.2 급료조사 데이터

```{r}
data_5.1<-read.table("All_Data/p130.txt",header=T,sep="\t")

# S:급료 X:경력 E:교육수준 M:관리(형태) / E,M은 범주형 변수

names(data_5.1)

# 범주형 질적 변수를 수치형으로 변형시켜서 예측하는데 사용한것이 질적 예측변수이다 

# E_1,E_2,E_3이런식으로 나누어서 0,1로 분류를 한다 (이것이 가변수)

# 더미변수를 만들 경우에는 역행렬의 조건에 의해서 -1개의 변수만 만들면 된다 

# 이는 공산성의 문제또한 있기에 이를 위해서 -1를 한것임 

### 자료형 변경 : 정수 -> 범주

data_5.1$E<-as.factor(data_5.1$E)

data_5.1$M<-as.factor(data_5.1$M)

head(data_5.1)

data_5.1$E #Levels: 1 2 3이라는 것이 생김 문자로 처리한다는 의미 

### 가변수 생성

data_5.1$E<-factor(as.character(data_5.1$E),levels = c("3","1","2")) 

#3번을 베이스 카테고리로 쓰기위한 설정 / 설정을 안하면 베이스는 e_1이 된다

data_5.1$M<-factor(as.character(data_5.1$M),levels = c("0","1")) 

data_dummy<-dummy_cols(data_5.1,

                       select_columns = c("E","M"),

                       remove_first_dummy = T,

                       remove_selected_columns = T) #첫번째 생성되는 더미 변수를 제거

data_dummy #가변수의 더미는 n-1개를 하는 것이 역행렬을 위한 것이기에 지워준다 

# 더미를 만든 이후 더미의 모체 변수인 E M을 지워주어야 한다 

### 회귀분석(1) - 가변수

res<-lm(S~.,data = data_dummy)

res 

### 회귀분석(2) - lm()

res_1<-lm(S~.,data=data_5.1)

res_1 

#더미변수를 많이 쓰기에 factor로 바꾸어주고 분석하면 알아서 더미변수를 만들어서 진행함

summary(res_1)

#E_3와 M_0는 Intercept에 포함이 되어있음 그렇기에 베이스 카테고리라고 한다 

#E가 3이고 M이 0이면 대학원이상 일반관리 직급 -> Intercept(Beta_0) + X*Beta_1 = Y

### 표준화잔차 대 경력연수

plot(data_5.1$X, rstandard(res_1),pch=19,xlab="범주",ylab="잔차")

# 0을 중심으로 잘 퍼져있는가를 확인해야함 

#### 표준화잔차 대 교육수준 - 관리 조합

EM<-paste0(data_5.1$E,data_5.1$M)

plot(EM,rstandard(res_1),pch=19,xlabl="범주",ylab="잔차")

### 상호작용 효과(Interaction Effect)

res<-lm(S~X+E+M+E*M,data=data_5.1)

res

summary(res)

### 표준화잔차 대 경력연수

plot(data_5.1$X,rstandard(res),pch=19,xlab="범주",ylab="잔차")

### Cook의 거리

plot(res,4) #33번째 데이터만 제외하고 다시 회귀모형을 생성예정

### 상호작용 효과 - 관측개체 33 제외 

data_use<-data_5.1[-33,]

res<-lm(S~X+E+M+E*M,data=data_use)

res

summary(res)

### 표준화잔차 대 경력연수

plot(data_use$X,rstandard(res),pch=19,xlab="범주",ylab="잔차")

#### 표준화잔차 대 교육수준 - 관리 조합

EM<-paste0(data_use$E,data_use$M)

plot(EM,rstandard(res),pch=19,xlabl="범주",ylab="잔차")

#상호 호과가 들어간 이 모형이 더 괜찮은 모양이라고 판단이 된다 

### 기본 급료 추정 - 표 5.6

res<-lm(S~X+E+M+E*M,data=data_use)

df_new<-data.frame(X=rep(0,6),

                   E=rep(1:3,c(2,2,2)),

                   M=rep(c(0,1),3))

### 가변수 생성 - 분석용

df_new$E<-factor(as.character(df_new$E),levels = c("3","1","2")) 

df_new$M<-factor(as.character(df_new$M),levels = c("0","1")) 

cbind(df_new,predict = predict(res,df_new,interval = "confidence"))

# 이를 보고 각 레벨에 따른 차이를 보고 얼마나 나는 지 분석이 가능해야 한다 

# ex) 고졸과 대학원졸의 관리자 직급의 급여의 차이는?(평균적으로)
```

### 5.4 회귀방정식의 체계: 두 집단의 비교

### 표 5.7 - 고용전 검사 프로그램 데이터

```{r}
data_5.7<-read.table("All_Data/p140.txt",header=T,sep="\t")

head(data_5.7)

# 모형 1 - 통합모형 인종간 차이가 없을때

model_1<-lm(JPERF~TEST,data=data_5.7)

summary(model_1) 

# 모형 3 

model_3<-lm(JPERF~TEST+RACE+TEST*RACE,data=data_5.7)

summary(model_3)

# 인종적인 차이가 있는지 없는지를 확인해야 한다 어떤 모형을 사용할지 

### H_0:gamma=delta=0

anova(model_1,model_3) #model_3가 FM(완전모형)

#P-value(0.05424)<0.05이기에 H_0 이기에 모형1을 선택하는 것이 옳다고 판단(그러나 확신x)

#서로의 R-squared롤 보면 model_3가 더 좋음 / ANOVA는 참고용 절대적이지는 않다 
```

### 그림 5.7 - 표준화잔차 대 검사점수 : 모형 1

```{r}
plot(data_5.7$TEST,rstandard(model_1),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.7 - 표준화잔차 대 검사점수 : 모형 1")
```

### 그림 5.8 - 표준화잔차 대 검사점수 : 모형 3

```{r}
plot(data_5.7$TEST,rstandard(model_3),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.8 - 표준화잔차 대 검사점수 : 모형 3")

# 통계에서 나오는 결과는 결정의 보조수단이지 절대적이지 않아 

# 3번 모형을 선택한다고 결정한다고 진행 

summary(model_3) # Multiple R-squared:  0.6643
```

### 그림 5.9 - 표준화잔차 대 검사점수 : 모형 1

```{r}
plot(data_5.7$RACE,rstandard(model_1),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.9 - 표준화잔차 대 검사점수 : 모형 1")

# 분리된 회귀분석 결과

data_5.7_R1<-subset(data_5.7,RACE==1)

model_R1<-lm(JPERF~TEST,data=data_5.7_R1)

summary(model_R1) #소수민족

data_5.7_R0<-subset(data_5.7,RACE==0)

model_R0<-lm(JPERF~TEST,data=data_5.7_R0)

summary(model_R0) #백인

# 통합모형에서 나온 각각의 회귀식이 통합모형에서 나온것과 동일함 따라서 인종별로 나누어서

# 진행할 필요없이 통일모형을 사용해서 진행을 하면 된다(이는 데이터셋을 나눈경우와 동일함)
```

### 그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만

```{r}
plot(data_5.7_R1$TEST,rstandard(model_R1),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.10 - 표준화잔차 대 검사점수 : 모형 1. 소수민족만만")
```

### 그림 5.11 - 표준화잔차 대 검사점수 : 모형 1 . 백인만

```{r}
plot(data_5.7_R0$TEST,rstandard(model_R0),

     pch=19,xlab="검사점수",ylab="잔차",

     main="그림 5.11 - 표준화잔차 대 검사점수 : 모형 1. 백인만")

# 적절한 합격점수의 결정 - 소수민족

# 고용전 검사점수의 합격점에 대한 95% 신뢰구간

ym<-4

xm<-(ym-0.09712)/3.31095

s<-1.292

n<-10

t<-qt(1-0.05/2,8)

c(xm-(t*s/n)/3.31095, xm+(t*s/n)/3.31095) #신뢰구간
```

### 5.4.2 동일한 기울기와 다른 절편항을 가지는 모형

```{r}
model_3<-lm(JPERF~TEST+RACE,data=data_5.7)

summary(model_3)

#Intercept = BETA_0 / TEST = BETA_1 / RACE = gamma

## H_0: gamma=0

anova(model_1,model_3) #p-value(0.1552)>0.05 이기에 H_0은 참이다//gamma=0

#기울기가 같고 절편이 다른 모형은 아니라고 데이터가 이야기하고 있다 

# 소수민족(RACE=1): (0.6120+1.0276)+2.2988*TEST = 1.6396+2.2988*TEST

# 백인(RACE=0): (0.6120)+2.2988*TEST
```

### 5.4.3 동일한 절편항과 다른 기울기를 가지는 모형

```{r}
model_3<-lm(JPERF~TEST+I(TEST*RACE),data=data_5.7)

summary(model_3) #I()를 하면 교호작용하는 것만 보이게 하려고 없으면 RACE항이 자동추가됨

## H_0: delta=0

anova(model_1,model_3) #p-value(0.03395)<0.05보다 작기에 delta항은 필요한 변수라는 사실 

#Intercept = BETA_0 / TEST = BETA_1 / I(TEST*RACE) = delta
```

## 6장 변수변환 : Transformation of Varibables

```{r}
library(dplyr)
```

### 6.3 x-선 방사에 의한 박테리아 사망률

```{r}
data_6.2<-read.table("All_Data/p168.txt",header=T,sep="\t")

head(data_6.2)
```

### 그림 6.5

```{r}
plot(N_t~t,data=data_6.2,pch=19)
```

### 6.3.1 선형모형의 부적절성

```{r}
res<-lm(N_t~t,data=data_6.2)

summary(res)
```

### 그림 6.6

```{r}
plot(data_6.2$t,rstandard(res),pch=19) # 표준화 잔차의 플랏

# 잔차가 랜덤하게 잘 퍼져있어야 하는데 적절한 회귀모형이 아니라는 사실이 나옴
```

### 6.3.2 선형성을 위한 로그변환

### 그림 6.7

```{r}
plot(log(N_t)~t,data=data_6.2,pch=19)

# 박테리아의 수에 로그를 취하니 선형성이 보인다 

res<-lm(log(N_t)~t,data=data_6.2)

summary(res)

plot(data_6.2$t,rstandard(res),pch=19)

# 로그를 취해주니 잔차가 랜덤하게 이루어져 있음 

# n_0에 대한 추론

exp(5.973160)

exp(coef(res)[1]) #로그를 취해주고 하는 부분이 잘 이해가 안됨 

exp(coef(res)[1]-0.0588/2) #불편추정량 구하기 (참고용)
```

### 6.4 분산안정화 변환

```{r}
data_6.6<-read.table("All_Data/p174.txt",header=T,sep="\t")

data_6.6 #운항률과 사고 발생 건수에 대한 자료

plot(Y~N,data=data_6.6,pch=19) #잔차의 분산이 계속 커지는 효과를 보임 

res_1<-lm(Y~N,data=data_6.6)

summary(res_1)

# 표준화잔차 대 N의 플롯 / 그림 6.11

plot(data_6.6$N,rstandard(res_1),pch=19) #등분산성에 만족하지 못하는 모형을 보인다 

# N에 대한 sqrt(Y)의 회귀

# N에 대한 Y의 회귀

res_2<-lm(sqrt(Y)~N,data=data_6.6)

summary(res_2)
```

### 그림 6.12

```{r}
plot(data_6.6$N,rstandard(res_2),pch=19) #0을 중심으로 잔차가 보다 잘 퍼져있음이 보임 
```

### 6.5 이분산성의 검출

```{r}
data_6.9<-read.table("All_Data/p176.txt",header=T,sep="\t")

data_6.9

# Y 대 X의 플로

plot(Y~X,data=data_6.9,pch=19,main="그림 6.13")

# X에 대한 Y의 회귀

res_1<-lm(Y~X,data=data_6.9)

summary(res_1)

# 표준화잔차 대 X의 플롯

plot(data_6.9$X, rstandard(res_1),pch=19,main = "그림 6.14")
```

### 6.6 이분산성의 제거

```{r}
# 변환된 Y/X와 1/X를 적합한 회귀

data_6.9_1<-data.frame(Y=data_6.9$Y/data_6.9$X,

                       X=1/data_6.9$X)

res_2<-lm(Y~X,data=data_6.9_1)

summary(res_2) #지금은 변환하고 프라임 값들의 추정치가 나온것이기에 원래 회귀식으로 돌아가야함 

# 본래 변환시 X를 나눴기에 X를 곱해준다 -> B_0,B_1의 추정값이 바뀐다 서로 

plot(data_6.9_1$X, rstandard(res_2),pch=19,

     xlab="1/X",ylab="잔차",main = "[그림 6.15]")

```

### 6.7 가중최소제곱법

```{r}
wt<-1/data_6.9$X^2 #가중값

res_3<-lm(Y~X,data=data_6.9,weights = wt)

summary(res_3) #결과는 위의 6.6과 동일한 값이 나옴 

```

### 6.8 데이터에 대한 로그변환

```{r}
data_6.9<-read.table("All_Data/p176.txt",header=T,sep="\t")

head(data_6.9)

# log(Y) 대 X의 산점도

plot((Y)~X,data=data_6.9,pch=19,main="식 6.9")

plot(log(Y)~X,data=data_6.9,pch=19,main="그림 6.16")

res_4<-lm(log(Y)~X,data=data_6.9)

summary(res_4)

# X에 대한 log(Y)의 회귀로 부터 얻은 표준화잔차플롯

plot(data_6.9$X,rstandard(res_4),pch=19,

     xlab="X",ylab="잔차",main="그림 6.17")

# X와 X^2에 대한 log(Y)의 회귀

df<-data.frame(log_Y=log(data_6.9$Y),

               X=data_6.9$X,

               X2=data_6.9$X^2)

res_5<-lm(log_Y~X+X2,data=df) #그러나 이러한 과정은 필요업고 아래방식으로..

res_5<-lm(log(Y)~X+I(X^2),data=data_6.9)

summary(res_5)

# 표준화잔차 대 적합값 플롯 : X와 X^2에 대한 log(Y)의 회귀

plot(fitted(res_5),rstandard(res_5),pch=19,

     xlab="X",ylab="잔차",main="그림 6.18")

# 표준화잔차 대 X의 플롯 : X와 X^2에 대한 log(Y)의 회귀

plot(data_6.9$X,rstandard(res_5),pch=19,

     xlab="X",ylab="잔차",main="그림 6.19")

# 표준화잔차 대 X^2의 플롯 : X와 X^2에 대한 log(Y)의 회귀

plot(data_6.9$X^2,rstandard(res_5),pch=19,

     xlab="X",ylab="잔차",main="그림 6.20") #2차항을 추가함으로서 더 잘 피팅됨 

# 7장은 스킵 / 8장 스킵
```

### 9.2 통계적 추론에 미치는 효과

```{r}
data_9.1<-read.table("All_Data/p236.txt",header=T,sep="\t")

head(data_9.1)

res<-lm(ACHV~.,data_9.1)

summary(res)

#F-statistic: 5.717 on 3 and 66 DF, p-value: 0.001535

#이는 의미있는 beta가 존재한다라는 의미

#그러나 각 계수들의 회귀계수를 보니 모두 0이라는 결과가 나옴 이는 다중공선성이다

# Correlation panel

panel.cor<-function(x,y){

par(usr=c(0,1,0,1))

r<-round(cor(x,y),digits = 3)

text(0.5,0.5,r,cex=1.5)

}

pairs(data_9.1[-1],lower.panel = panel.cor)

```

### 그림 9.1

```{r}
plot(fitted(res),rstandard(res),pch=19,

xlab="예측값",ylab="잔차",main="[그림 9.1]")

#이는 회귀모형에 동시에 들어가면 안된다라는 의미 (제거가 필요)
```

### 9.3 예측에 미치는 효과

```{r}
data_9.5<-read.table("All_Data/p241.txt",header=T,sep="\t")

head(data_9.5)

### 산점도

pairs(data_9.5[-1],lower.panel = panel.cor)

# 회귀분석(1) : 데이터 1949~1966

res<-lm(IMPORT~.,data_9.5[-1])

summary(res) #다중공선성이 존재한다고 예측을 일단함
```

### 그림 9.3 : 표준화잔차의 인덱스플롯

```{r}
plot(1:nrow(data_9.5),rstandard(res),

pch=19,type="b",

xla="번호",ylab="잔차",main="[그림 9.3]")

# 회귀분석(2) : 데이터 1949~1959

data_use<-subset(data_9.5,YEAR<=59)

res<-lm(IMPORT~.,data_use[-1])

summary(res)

```

### 그림 9.4 : 잔차플롯

```{r}
plot(1:nrow(data_use),rstandard(res),

pch=19,type="b",

xla="번호",ylab="잔차",main="[그림 9.4]")
```

### 9.4 다중공선성의 탐색

```{r}
# 분산확대인자

library(olsrr)

# 교육기회 균등(EEO)

res_9.1<-lm(ACHV~.,data_9.1)

summary(res_9.1)

#car::vif(res_9.1)

olsrr::ols_vif_tol(res_9.1) # 10보다 크면 유의한 영향을 준다 제거필요?

```

### 표 9.5 : 프랑스 경제 데이터

```{r}
data_use<-subset(data_9.5,YEAR<=59)

res_9.5<-lm(IMPORT~.,data_use[-1])

#car::vif(res_9.1)

olsrr::ols_vif_tol(res_9.5) #작으면 새로운 변수를 만들거나 제거를 통해서 진행

# 상태 지수

cnd.idx<-olsrr::ols_eigen_cindex(res_9.1)

round(cnd.idx,4)

cnd.idx<-olsrr::ols_eigen_cindex(res_9.5)

round(cnd.idx,4)
```

## 10장 - 공선형 데이터의 처리

```{r}
data_9.5<-read.table("All_Data/p241.txt",header=T,sep="\t")

head(data_9.5)

 

# 회귀분석(2) : 데이터 1949~1959

data_use<-subset(data_9.5,YEAR<=59)

res<-lm(IMPORT~.,data_use[-1]) #1열 제거 

summary(res)

head(data_use[-c(1:2)])

### 주성분(principle component)

pc<-prcomp(data_use[-c(1:2)],scale.=T)

pc$rotation

pc$x #변화된 새로운 변수가 저장된곳

### 주성분회귀 - 원래데이터를 주성분분석을 통해 새로운 데이터를 생성 이를 가지고 회귀 

df<-data.frame(IMPORT = scale(data_use$IMPORT),

               pc$x)

df

res<-lm(IMPORT~.,data=df)

summary(res) #유의한 1,2번째 계수들만 사용하고 3번째 것은 없이 해도 되겠다..

#중간고사 이후의것이 주가나오겠지만 앞에것도 알아야함 연습문제에서 나올것 예상

```

## 11장

```{r}
# 11.10 - 

data_3.3<-read.table("All_Data/p060.txt",header=T,sep="\t")

head(data_3.3)

# 분산확대 인자(VIF) : 10초과 > 심각한 공산성의 문제가 있음

res<-lm(Y~.,data=data_3.3)

summary(res)

olsrr::ols_vif_tol(res) #10을 넘는 것이 없기에 공산성의 문제는 없다 

# 수정결정계수

res_summ<-summary(res)

res_summ$adj.r.squared #0.662846

# Mallow's Cp - 작을수록 좋음

# 축소모형

res_subset<-lm(Y~X1+X3,data=data_3.3)

olsrr::ols_mallows_cp(res_subset,res) 

# AIC - 참고만

olsrr::ols_aic(res_subset,method = "SAS")

# BIC

olsrr::ols_sbic(res_subset,res)

### 전진적 선택법 - AIC

res_step<-step(res,direction = "forward")

### 후진적 제거법 -AIC

res_step<-step(res,direction = "backward")

### 단계적 방법

res_step<-step(res) #최종적으로는 X1+X3이 선택됨 이거 사용하는 것이 좋을듯 

summary(res_step) #X3가 의미 없다고 나와도 아님 검증된것임
```

## 0613 - regression analysis

### 연습문제 11.5

```{r}
data_use<-read.table("All_Data/p329.txt",header=T,sep="\t")

head(data_use)

### 데이터 탐색 - 자료형

sapply(data_use,class)

### 데이터 탐색 - 기초 통계량

summary(data_use) #결측값의 여부를 확인

### 데이터 탐색 - 히스토그램 & 상자그림

hist(data_use$Y,main="histogram of data_use$Y")

boxplot(data_use$Y)

### 데이터 탐색 - 산점도 행렬 & 상관계수

plot(data_use)

cor(data_use) 

pairs(data_use)

panel.cor<-function(x,y){

  par(usr=c(0,1,0,1))

  r<-round(cor(x,y),digits = 3)

  text(0.5,0.5,r,cex=1.5)

}

pairs(data_use,lower.panel = panel.cor)

### (a) 모든 변수들이 모형에 포함시킬 것인가?

cor(data_use)

model_1<-lm(data_use$Y~.,data_use)

summary(model_1) #전형적인 다중공선성이 보인다 

#coef가 X6 0.820이므로 확인 필요

olsrr::ols_vif_tol(model_1) #X6,X7이 10에 가깝기에 제거해야하는 대상중 최우선 

#다중공선성이 보이기에 모든 변수를 모형에 포함 시킬수는 없다 

#10보다 크면 심각한 공선성이 존재 

### (b) 지방세(X1) 방의 수(X6), 건물의 나이(X8)가 판매가격(Y)을 

#       설명하는데 적절하다는 의견에 동의 하는가?

model_2<-lm(Y~X1+X6+X8,data=data_use)

summary(model_2)

# VIF

olsrr::ols_vif_tol(model_2)

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_2)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_2) 

#동의는 할 수 있으나(적절하지만) 최고의 모형이라는데는 동의 못함

### (c) 지방세 X1가 단독으로 판매가격 Y을 설명하는데 적절하다는 의견에 동의?

model_3<-lm(Y~X1,data=data_use)

summary(model_3) #adj-rsquared는 변수를 선택할때 사용 

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_3)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_3) 

### 적절한 모형 제시

data_use_2<-data_use[-6]

model_4<-lm(Y~.,data_use_2)

summary(model_4)

# VIF

olsrr::ols_vif_tol(model_4)

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_4)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_4) 

### 단계적 선택법 - AIC

model_5<-step(model_4)

summary(model_5)

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_5)

layout(1)

# Cook의 거리 - 0.5보다 작으면 괜찮음 / 1보다 큰 값을 고려 / 전체적으로 봤을때 튀는 값이 있는 경우 

olsrr::ols_plot_cooksd_chart(model_5) 

### 17번 제거

data_use_3<-model_5$model[-17,]

model_6<-lm(Y~.,data_use_3)

summary(model_6)

# VIF

olsrr::ols_vif_tol(model_6)

# 회귀진단 그래프들

layout(matrix(1:4,nrow=2,byrow=T))

plot(model_6)

layout(1)

# Cook의 거리 - 1을 넘으면 확인을 해야한다 

olsrr::ols_plot_cooksd_chart(model_6)
```